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

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

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.

_, 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()

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

from scipy.io.wavfile import write
 
# SAMPLE_RATE = 44100
write("sinewave.wav", SAMPLE_RATE, mix_norm)

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

# 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()

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

from scipy.fft import irfft
 
new_sig = irfft(yf)
 
plt.plot(new_sig[:1000])
plt.show()