====== Signal processing ======

====== 1. Harmonický signál ======

V první ukázce vygenerujeme několik period harmonického signálu o frekvenci 1000 Hz. Pokud bychom chtěli data pohodlně přehrát na počítači, uložíme je do souboru typu ''wav''. 

Asi nikoho nepřekvapí matematický popis časového průběhu harmonického signálu s frekvencí $f$ a amplitudou $A$ 

$y(t) = A * sin(2 * \pi * f * t)$

Protože nás zajímá digitální (tj. vzorkovaná) podoba signálu, je třeba se zabývat vlastnostmi vzorků:

  * vzorkovací frekcence musí odpovídat [[https://cs.wikipedia.org/wiki/Nyquist%C5%AFv%E2%80%93Shannon%C5%AFv_vzorkovac%C3%AD_teor%C3%A9m|Nyquistově podmnínce]]. S ohledem na maximální slyšitelnou frekcenci cca 20 kHz můžeme zvolit vzorkovací frekvenci 44100 Hz, která je obvyklá v audiotechnice
  * amplituda v případě uložení do souboru nemůže být být neceločíselného datového typu ''float''. Datový formát ''wav'' podporuje pouze celá kladná čísla v rozsahu 16 bit;, takže je třeba provést přeškálování (normalizaci). 

<code python>
import numpy as np
from matplotlib import pyplot as plt

SAMPLE_RATE = 44100  # vzorkovaci frekvence v Hz
DURATION = 5         # doba trvani

def generate_sine_wave(freq, sample_rate, duration):
    x = np.linspace(0, duration, sample_rate * duration, endpoint=False)
    frequencies = x * freq
    # argumentem np.sin jsou radiany 
    y = np.sin((2 * np.pi) * frequencies)
    return x, y

# 5s signalu 1kHz 
x, y = generate_sine_wave(2, SAMPLE_RATE, DURATION)
plt.plot(x, y)
plt.title('Harmonicky signal 2Hz')
plt.xlabel('Cas (s)')
plt.ylabel('Amplituda')
plt.show()
</code>

===== Vytvoření zašuměného signálu =====

Pro další experiment budeme simulovat harmonický signál obsahující šum: vytvořme dva různé signály se stejnou časovou osou a rozdílnými frekvencemi, amplitudu šumového signálu potlačíme. Abychom mohli signál uložit do souboru ''wav'', provedeme normalizaci amplitudy.

<code python>
_, signal = generate_sine_wave(400, SAMPLE_RATE, DURATION)
_, sum = generate_sine_wave(4000, SAMPLE_RATE, DURATION)
sum = sum * 0.3

mix = signal + sum

mix_norm = np.int16((mix / mix.max()) * 32767)

plt.plot(mix_norm[:1000])
plt.show()
</code> 

K uložení do souboru použijeme vestavěnou podporu knihovny scipy

<code python>
from scipy.io.wavfile import write

# SAMPLE_RATE = 44100
write("sinewave.wav", SAMPLE_RATE, mix_norm)
</code>

===== Analýza pomocí Fourierovy transformace =====

<code python>
from scipy.fft import fft, fftfreq

# Number of samples in normalized_tone
N = SAMPLE_RATE * DURATION

yf = fft(mix_norm)
xf = fftfreq(N, 1 / SAMPLE_RATE)

plt.plot(xf, np.abs(yf))
plt.show()
</code>

===== Filtrace signálu ve spektru =====

<code python>
# maximální frekvence je právě na polovině vzorkovací fr.
points_per_freq = len(xf) / (SAMPLE_RATE / 2)

# index v seznamu odpovídající 4000 Hz
target_idx = int(points_per_freq * 4000)

# vynulování všech položek seznamu, kromě frekvence 4k
yf[target_idx - 1 : target_idx + 2] = 0

plt.plot(xf, np.abs(yf))
plt.show()
</code>

Inverzní Fourierovou transformací pak můžeme získat filtrovaný signál

<code python>
from scipy.fft import irfft

new_sig = irfft(yf)

plt.plot(new_sig[:1000])
plt.show()
</code>

<note tip>Popisovaná metoda není příliš vhodná na reálné signály, protože v praxi nebudeme mít téměř nikdy čistý harmonický signál a je třeba pouřít vhodný filtr.</note>

====== 2. Filtrace signálu pomocí FIR filtru ======