<jupyter_start><jupyter_text>Fourierova integralska transformacija Fourierova transformacija je matematično orodje, s katerim preslikamo signale **iz časovne v frekvenčno domeno**. :::{note}Fourierova transformacija in njen inverz sta integralski preslikavi, ki ste ju definirali na predavanjih:$$X(f)=\int_{-\infty}^{+\infty} x(t)\,e^{-\textrm{i}\,2\pi\,f\,t}\,\textrm{d} t,\qquadx(t)=\int_{-\infty}^{+\infty} X(f)\,e^{\textrm{i}\,2\pi\,f\,t}\,\textrm{d} f.$$::: :::{note}Realni (pomerjeni) signali običajno niso podani v obiliki funkcijskega predpisa $x(t)$, temveč so podani v obliki časovne vrste vrdnosti signla, pomerjene pri diskretnih časovnih točkaj, $x(n\,\Delta t)$.V praksi bomo zato pogosto raje govorii o **Diskretni Fourierovi transformaciji (DFT)**, definirani na predpostavki diskretno vzorčenega s konstantnim časovnim korakom, katere rezultat je *diskretna frekvenčna prdstavitev signala s konstantnim korakom v frekvenčni domeni*, $X(k\,\Delta f)$. Na področju procesiranja signalov je izrednega pomena tudi numerično učinkovit algoritem **Hitre Fourierove transformacije (FFT)**, ki omogoča časovno-frekvenčno preslikavo zajetih diskretnih signalov v realnem času.O DFT in FFT boste več govorili na enem od naslednjih predavanj.:::  Fourierova transformacija v Pythonu Poglejmo si uporabo algortima (diskretne) Fourierove transformacije, implementiranega v Numpy in Scipy, na primeru sinusnega signala z amplitudo $A$ ter sekundno frekvenco $p$. Na predavanjih ste povedali, da velja:$$x(t) = A\,\sin(2\pi\,p\,t) = \frac{A}{2\,\textrm{i}}\,\left(e^{\textrm{i}\,2\pi\,p\,t)}-e^{-\textrm{i}\,2\pi\,p\,t)}\right)$$Fourierova transformacija $x(t)$ pa je:$$\mathcal{F}_{t}\left( A\,\sin(2\pi\,p\,t)\right)=\frac{\mathrm{i} A}{2} \delta(f - p) - \frac{\mathrm{i} A}{2}\delta(f + p)$$ Pripravimo funkcijo za vizualizacijo Fourierove transformacije sinusnega signala:<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
def fourier_sin(A, p):
    """
    Vizualizacija fourierove transformacije sinusnega signala z amplitudo
    `A` in frekvenco `p` Hz.
    
    Vrne amplitudo spektra pri frekvenci `p`.
    """
    A_p = 1j*A/2
    
    freq = np.array([-2*p, -p, 0, p, 2*p])
    amp = np.array([0., A_p, 0, A_p, 0])
    
    plt.axhline(y=0, c='k', lw=2)
    plt.axvline(x=0, c='k', lw=1)

    for i, x in (enumerate(freq)):
        plt.plot([x, x], [0, np.abs(amp[i])], 'k', lw=2)
    plt.plot(freq, np.abs(amp), 'ro', label='analitični amplitudni spekter')
    
    plt.xlabel('$f [Hz]$')
    plt.ylabel('$|X(f)|$')

p = 5
A = 10
fourier_sin(A, p)<jupyter_output><empty_output><jupyter_text>Pripravimo sinusni signal in poglejmo rezultat DFT:<jupyter_code>dt = 1/(p*50)
t = np.arange(0, 1, dt)
x = A*np.sin(2*np.pi*t*p)<jupyter_output><empty_output><jupyter_text>Izračun FFT poljubnega signala:<jupyter_code>X = np.fft.fft(x)
freq = np.fft.fftfreq(len(t), dt)<jupyter_output><empty_output><jupyter_text>Izračun FFT realnega signala:<jupyter_code>X_r = np.fft.rfft(x)
freq_r = np.fft.rfftfreq(len(t), dt)
fourier_sin(A, p)
plt.plot(freq, np.abs(X), '.-', label='Numpy FFT poljubnega signala')
plt.plot(freq_r, np.abs(X_r), '.-', label='Numpy FFT realnega signala')
plt.xlim(-2*p, 2*p)
plt.legend();<jupyter_output><empty_output><jupyter_text>Opzimo, da se teoretične in dobljene amplitude ne ujemajo. :::{admonition} Naloga 1 (10 minut):class: importantRaziščite implementacijo skaliranja (normalizacije) [diskretne Fourierove transformacije v Numpy](https://numpy.org/doc/stable/reference/routines.fft.html) in zgornjo kodo nadgradite tako, da se bodo dobljene amplitude ujemale s teoretičnimi.*Dodatek*: Raziščite možnosti povečanja frekvenčne ločljivosti ($\downarrow \Delta f$) pri uporabi diskretne fourierove transformacije. S katerim parametrom zajema signala je ta poezana? (Namig: argmument `n` funkcije `np.fft.rfft`.):::<jupyter_code>ločljivost = 10
X_1 = np.fft.rfft(ločljivost*x, norm='forward', n=len(t)*ločljivost)
freq_ = np.fft.rfftfreq(n=len(t)*ločljivost, d=dt)
fourier_sin(A, p)
plt.plot(freq_, np.abs(X_1), '.-', label='Numpy FFT realnega signala')
plt.xlim(-2*p, 2*p)
plt.legend()
plt.title('Možnost 1: frekvenčna interpolacija');
t_2 = np.arange(0, ločljivost, dt)
x_2 = A*np.sin(2*np.pi*t_2*p)
X_2 = np.fft.rfft(x_2, norm='forward', n=len(t_2))
fourier_sin(A, p)
plt.plot(freq_, np.abs(X_2), '.-', label='Numpy FFT realnega signala')
plt.xlim(-2*p, 2*p)
plt.legend()
plt.title('Možnost 2: sprememba parametrov zajema');<jupyter_output><empty_output><jupyter_text>DFT realnega signala (03-naloga-2)=:::{admonition} Naloga 2 (30 min):class: importantZ uporabo generatorja signalov in Arduino zajemnega sistema pripravite ter zajamite dva signala različnih oblik s parametri iz podatkov domače naloge 3, ki jih najdete v [e-učilnici](https://e-ucilnica.fs.uni-lj.si/course/view.php?id=328).Pri zajemu podatkov lahko uporabite LabView program, ki ste ga pripravili na prejšnji vaji, ali pa program, ki ga najdete spodaj.Ne pozabite programa pred zajemom prirediti tako, da bo:* deloval z vašim zajemnim sitemom (nastavitev COM vrat, mesto shranjevanja .lvm datoteke),* če program omogoča nastavitev občutljivosti, te nastavite na vrednost 1,* omogočal nastavitev časa trajana segmenta (namesto števila vzorcev v segmentu),* zajel signal glede na parametre naloge (pozor, tokrat shranjujemo le en segment meritve).Nalogo lahko rešujete v parih, tako, da si ustrezno razdelite delo (nastavitev generatorja signala, zajem signalov).Osnovna navodila za uporabo generatorja signalov RIGOL DG 1022 najdete v [predlogi prejšnje laboratorijske vaje](02-uporaba-rigol).::: ![labview-png](./images/03/zajem_N_channels_closed_loop.png)Za delovanje potrebujete tudi naslednji podprogram:![labview-podprogram](./images/03/podprogram_clip_fs.png)Oboje lahko prenesete lahko tudi v obliki {download}`zip arhiva`.  Obdelava zajetega signala :::{admonition} Naloga 3 (15 min):class: seealso* Sinusni signal, zajet pri prejšnji nalogi, naložimo v Python program.* Določimo parametre (amplitudo $A$, frekvenco $p$, fazni zamik $\varphi$) zajetega signala (`scipy.optimize.curve_fit`).* Amplitudni spekter zajetega sinusnega signala primerjajmo s teoretičnim amplitudnim spektrom funkcije $A \, \sin(2\, \pi \, p \, t)$.::: Pri branju podatkov, shranjnih v `.lvm` datoteke, si lahko pomagate s paketom [`lvm_read`](https://pypi.org/project/lvm-read/).Primer uporabe si poglejmo spodaj.<jupyter_code>import lvm_read
file = 'data/03/sinus.lvm'
data = lvm_read.read(file)
data['Segments']<jupyter_output><empty_output><jupyter_text>Segment podatkov je v slovarju shranjen s ključem zaporednega indeksa. Vsak posameni segment je prav tako slovar, v katerem do surovih shranjenih podatkov dostopamo s ključem `'data'`:<jupyter_code>kanali = data[0]['data']
kanali.shape<jupyter_output><empty_output><jupyter_text>Posamezni stolpci podatkov predstavljajo po en shranjen kanal:<jupyter_code>signal = kanali[:, 0]
signal.shape<jupyter_output><empty_output><jupyter_text>Časovni vektor lahko pripravimo na podlagi znane frekvence vzorčenja:<jupyter_code>fs = 100 # vzorcev / s
t = np.arange(len(signal)) / fs
plt.figure()
plt.plot(t, signal, label='kanal 1')
plt.xlabel('t [s]')
plt.ylabel('napetost [V]')
plt.legend()<jupyter_output><empty_output><jupyter_text>Določitev parametrov zajetega signala s `scipy.optimize.curve_fit`:<jupyter_code>from scipy.optimize import curve_fit<jupyter_output><empty_output><jupyter_text>Definiramo funkcijo modela našega signala, katere prvi parameter je neodvisna spremenljivka $t$, ostali argumenti pa so iskani parametri:<jupyter_code>def model(t, A, p, phi, V_0):
    """
    Model funkcije `A*sin(2*pi*p*t + phi) + V_0` za aproksimacijo.
    """
    return A*np.sin(2*np.pi*p*t + phi) + V_0<jupyter_output><empty_output><jupyter_text>Določimo smiselne začetne približke:<jupyter_code>A_0 = (np.max(signal) - np.min(signal)) / 2
plt.plot(t, signal, 'k')
plt.fill_between(t, np.mean(signal) - A_0, np.mean(signal) + A_0, color='C0', alpha=0.25)
plt.xlabel('t [s]')
plt.ylabel('signal [V]');
from scipy.signal import find_peaks
i_vrhov = find_peaks(signal)[0]
t_vrhov = t[i_vrhov]
p_0 = 1 / (t_vrhov[1] - t_vrhov[0])
plt.plot(t, signal, 'k')
plt.fill_between(np.linspace(t[i_vrhov[0]], t[i_vrhov[0]] + 1/p_0), np.min(signal), np.max(signal), color='C0', alpha=0.25)
plt.xlabel('t [s]')
plt.ylabel('signal [V]');
V_0 = np.mean(signal)
plt.plot(t, signal, 'k')
plt.axhline(V_0, c='r')
plt.xlabel('t [s]')
plt.ylabel('signal [V]');
normal_signal = (signal - V_0) / A_0
phi_0 = np.arcsin(normal_signal[0])
plt.plot(t, normal_signal, 'k')
plt.plot(t[0], normal_signal[0], 'ro')
print(np.rad2deg(phi_0))
plt.xlabel('t [s]')
plt.ylabel('normaiziran signal [V]');
popt, pcov = curve_fit(model, t, signal, p0=[A_0, p_0, phi_0, V_0])
A, p, phi, V_0 = popt
plt.plot(t, model(t, A, p, phi, V_0), lw=2, label='aproksimiran signal')
plt.plot(t, signal, 'k.', label='zajet signal')
plt.legend()
plt.xlabel('t [s]')
plt.ylabel('signal [V]');<jupyter_output><empty_output><jupyter_text>(Diskretna) Fourierova transformacija zajetega signala<jupyter_code>S = np.fft.fft(signal) / len(t)
freq = np.fft.fftfreq(len(t), 1/fs)
fourier_sin(A, p)
plt.plot(freq, np.abs(S), '.-', label='zajeti signal')
plt.xlim(-2*p, 2*p);
plt.legend(loc=(1.01, 0));<jupyter_output><empty_output><jupyter_text>Odprava statične (DC) komponte:<jupyter_code>S_1 = np.fft.fft(signal - V_0) / len(t)
freq = np.fft.fftfreq(len(t), 1/fs)
fourier_sin(A, p)
plt.plot(freq, np.abs(S_1), '.-', label='zajeti signal brez statične komponente')
plt.xlim(-2*p, 2*p);
plt.legend(loc=(1.01, 0));<jupyter_output><empty_output><jupyter_text>Pogosto pri realnih signalih izračunamo amplitudni spekter le pri pozitivnih frekvencah.Da se energija signala ohrani, dobljen enostranski spekter ustrezno skaliramo:<jupyter_code>S_r = np.fft.rfft(signal - V_0) / len(t)
S_r[1:] *= 2
freq = np.fft.rfftfreq(len(t), 1/fs)
fourier_sin(A, p)
plt.plot(freq, np.abs(S_r), '.-', label='enostranski spekter realnega signala')
plt.axhline(A, c='g', label='amplituda zajetega signala')
plt.xlim(-2*p, 2*p)
plt.legend(loc=(1.01, 0));
np.isclose(np.sum(np.abs(S_1)),  np.sum(np.abs(S_r)), rtol=1e-4)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15430586852160283, 'fraction_numerical': 0.01596943274223572, 'mean_word_length': 6.04, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 61, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '600221', 'license': 'MIT', 'n_tokens_mistral': 4574, 'n_tokens_neox': 4437, 'n_words': 934, 'path': 'lv_book/notebooks/03 - Fourierova integralska transformacija.ipynb', 'repo_name': 'domengorjup/ps_lv', 'size': '10207'}