Search
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.
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ů:
float
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()