Fourierova integralska transformacija

Fourierova integralska transformacija#

Domača naloga#

Domača naloga

  • Oba signala, zajeta pri Nalogi 2, odprite v izbranem programskem okolju (npr. Python).

  • Prikažite in primerjajte njuna amplitudna ter fazna spektra.

  • Enaka signala generirajte tudi numerično (na primer z uporabo scipy.signal). Primerjajte amplitudna in fazna spektra zajetih ter numerično pripravljenih signalov.

Pripravite kratko poročilo (obseg od 3 do 10 celic s kodo v okolju Jupyter Notebook), iz katerega naj bodo razvidni podatki naloge (iz tabele), ter da ste vse parametre in zahteve naloge naloge tudi upoštevali. Poročilo oddajte v .pdf obliki (glejte navodila za oddajo domačih nalog).

Dodatek: Raziščite lastnost časovnega premika Fourierove transformacije, in:

  • ocenite časovni zamik med zajetima in generiranima signaloma.

  • Zajeta signala v frekvenčni domeni z uporabo lastnosti časovnega zamika poravnajte z generiranima.

  • Poravnana signala preslikajte nazaj v časovno domeno (np.fft.irfft) in rezultat poravnave prikažite.

Note

Lastnost časovnega premika je opisana v predlogi predavanja 3. Časovni premik signala za vrednost \(t_0\) je v frekvenčni domeni opisan z:

\[ \mathcal{F}\{x(t-t_0)\} = e^{-\textrm{i}\,2\pi\,f\,t_0}\,X(f). \]

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, \qquad x(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 obliki funkcijskega predpisa \(x(t)\), temveč so podani v obliki časovne vrste vrednosti signala, pomerjene pri diskretnih časovnih točkah, \(x(n\,\Delta t)\).

V praksi bomo zato pogosto raje govorili o Diskretni Fourierovi transformaciji (DFT), definirani na predpostavki diskretno vzorčenega signala s konstantnim časovnim korakom, katere rezultat je diskretna frekvenčna predstavitev 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 algoritma (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:

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)
../_images/4d7ad14c49c7a48d23916c79fc72a07de97c4e28296b3a38f4fa4956f9672c51.png

Pripravimo sinusni signal in poglejmo rezultat DFT:

dt = 1/(p*50)
t = np.arange(0, 1, dt)
x = A*np.sin(2*np.pi*t*p)

Izračun FFT poljubnega signala:

X = np.fft.fft(x)
freq = np.fft.fftfreq(len(t), dt)

Izračun FFT realnega signala (\(x(t) \in \mathcal{R}\)):

X_r = np.fft.rfft(x)
freq_r = np.fft.rfftfreq(len(t), dt)

Hide code cell source

 
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();
../_images/f114fc8e09151ae6d35c3e683b4fac097251b98e35f67eb54b419b822d7c8eda.png

Opazimo, da se teoretične in dobljene amplitude ne ujemajo.

Naloga 1 (10 minut)

Raziščite implementacijo skaliranja (normalizacije) diskretne Fourierove transformacije v Numpy 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 povezana? (Namig: argument n funkcije np.fft.rfft.)

../_images/5ef5a432b5c363a5f43a586ae224304c61e019258250e71d47c2c1d722e411c4.png
../_images/c9dd06e3cc4d261101e46c2fca757de6e9c942bb98e45d7c306c3a1d02774c2a.png

DFT realnega signala#

Naloga 2 (30 min)

Z uporabo generatorja signalov in Arduino zajemnega sistema pripravite ter zajemite dva signala različnih oblik s parametri iz podatkov domače naloge 3, ki jih najdete v e-učilnici.

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 sistemom (nastavitev COM vrat, mesto shranjevanja .lvm datoteke),

  • če program omogoča nastavitev občutljivosti, te nastavite na vrednost 1,

  • omogočal nastavitev časa trajanja segmenta (namesto števila vzorcev v segmentu),

  • zajel signal glede na parametre naloge.

Nalogo lahko rešujete v grupah, 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.

labview-png

Za delovanje potrebujete tudi naslednji podprogram:

labview-podprogram

Oboje lahko prenesete lahko tudi v obliki zip arhiva.

Obdelava zajetega signala#

Naloga 3 (15 min)

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

Primer uporabe si poglejmo spodaj.

import lvm_read
file = 'data/03/sinus.lvm'
data = lvm_read.read(file)
data['Segments']

1

Segment podatkov je v slovarju shranjen s ključem zaporednega indeksa. Vsak posamezni segment je prav tako slovar, v katerem do surovih shranjenih podatkov dostopamo s ključem 'data':

kanali = data[0]['data']
kanali.shape

(500, 2)

Posamezni stolpci podatkov predstavljajo po en shranjen kanal:

signal = kanali[:, 0]
signal.shape

(500,)

Časovni vektor lahko pripravimo na podlagi znane frekvence vzorčenja:

fs = 100 # vzorcev / s
t = np.arange(len(signal)) / fs

Hide code cell source

 
plt.figure()
plt.plot(t, signal, label='kanal 1')
plt.xlabel('t [s]')
plt.ylabel('napetost [V]')
plt.legend()

<matplotlib.legend.Legend at 0x7f1d8a0b3c50>
../_images/2f60330f0c0ed86016ad8772cb093aeb6cf356f34219e7cb5487cb083fd02675.png

Določitev parametrov zajetega signala s scipy.optimize.curve_fit:#

from scipy.optimize import curve_fit

Definiramo funkcijo modela našega signala, katere prvi parameter je neodvisna spremenljivka \(t\), ostali argumenti pa so iskani parametri:

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

Določimo smiselne začetne približke:

A_0 = (np.max(signal) - np.min(signal)) / 2

Hide code cell source

 
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]');
../_images/99779254359508acd944de838c25dbec47c7ffa466287eedc0521618b5239c33.png 
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])

Hide code cell source

 
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]');
../_images/763a74ff1fdd3d8df6a40bb5542a2ad22acd1bc28dbc94527b01dc6965fb811d.png 
V_0 = np.mean(signal)

Hide code cell source

 
plt.plot(t, signal, 'k')
plt.axhline(V_0, c='r')
plt.xlabel('t [s]')
plt.ylabel('signal [V]');
../_images/1e60752266c203a44e22529e39e1613abbe936730d335ceba7357c1c1b164077.png 
normal_signal = (signal - V_0) / A_0
phi_0 = np.arcsin(normal_signal[0])

Hide code cell source

 
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('normaliziran signal [V]');

45.003549449490734
../_images/599a72a44d0e1baaa0ca02c43e5e8a7ab02d5ef76fcdf8352a888a4d3aa98ee6.png 
popt, pcov = curve_fit(model, t, signal, p0=[A_0, p_0, phi_0, V_0])
A, p, phi, V_0 = popt

Hide code cell source

 
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]');
../_images/6cf1cdb3f8159c903c18eca92a00bb6744c7e61bb77d2ebe7f6a699196259130.png

(Diskretna) Fourierova transformacija zajetega signala#

S = np.fft.fft(signal) / len(t)
freq = np.fft.fftfreq(len(t), 1/fs)

Hide code cell source

 
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));
../_images/19184b35c57d0a5032359eaff18b796f1fdae49090ce86ecf455f7ef993640ff.png

Odprava statične (DC) komponte:

S_1 = np.fft.fft(signal - V_0) / len(t)
freq = np.fft.fftfreq(len(t), 1/fs)

Hide code cell source

 
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));
../_images/e6e166b44ac9f1714b317028c2a222fe3c9bfc184fcf594af08881650425a89a.png

Pogosto pri realnih signalih izračunamo amplitudni spekter le pri pozitivnih frekvencah.

Da se energija signala ohrani, dobljen enostranski spekter ustrezno skaliramo:

S_r = np.fft.rfft(signal - V_0) / len(t)
S_r[1:] *= 2
freq = np.fft.rfftfreq(len(t), 1/fs)

Hide code cell source

 
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));
../_images/188182759541aa3c9499f17c8f7c49d6206be3a8af76b201820b6bcc02e0a410.png 
np.isclose(np.sum(np.abs(S_1)),  np.sum(np.abs(S_r)), rtol=1e-4)

np.True_

Amplitudni in fazni spekter realnega signala#

Poglejmo primer amplitudnega in faznega spektra realnega signala z dvema harmonskima komponentama.

Parametri signala:

A1 = 1
A2 = 2
f1 = 1
f2 = 2
phi1 = np.pi/2
phi2 = np.pi/4

Parametri vzorčenja:

fs = 50
dt = 1/fs
T = 5
t = np.arange(0, T*fs) * dt

Obravnavan signal:

s1 = A1*np.cos(2*np.pi*f1*t + phi1)
s2 = A2*np.cos(2*np.pi*f2*t + phi2)
signal = s1 + s2

Hide code cell source

 
plt.figure()
plt.plot(t, signal, label='signal = s1 + s2')
plt.plot(t, s1, label='s1')
plt.plot(t, s2, label='s2')
plt.xlabel('t [s]')
plt.ylabel('signal [/]')
plt.legend(loc=(1.01, 0));
../_images/e5e246fb19ecf5d18cb2c3b79cbf46defd24700d01a88f5c8b4f3e4dc47aedc1.png

Fourierova transformacija signala, normiranje enostranskega spektra:

S = np.fft.rfft(signal, norm='forward')
S[1:] *= 2
freq = np.fft.rfftfreq(len(t), dt)

Osnovne frekvence komponent \(s_1\) in \(s_2\) obravnavanega signala:

i1 = np.argmin(np.abs(freq-1))
i2 = np.argmin(np.abs(freq-2))
freq[i1], freq[i2]

(np.float64(1.0), np.float64(2.0))

Hide code cell source

 
fig, ax = plt.subplots(2, 1)
ax[0].plot(freq, np.abs(S), 'C0', label='amplitudni spekter')
ax[0].set_xlim(0, 5)
ax[0].set_xlabel('f [Hz]')
ax[0].set_ylabel('$|S|$ [/]')
ax[0].axvline(f1, c='k', lw=0.5)
ax[0].axvline(f2, c='k', lw=0.5)
ax[0].axhline(A1, c='k', lw=0.5)
ax[0].axhline(A2, c='k', lw=0.5)
ax[0].legend(loc=(1.01, 0))

ax[1].plot(freq, np.angle(S), 'C1', label='fazni spekter')
ax[1].set_xlim(0, 5)
ax[1].set_xlabel('f [Hz]')
ax[1].set_ylabel('$\\angle S$ [rad]')
ax[1].axvline(f1, c='k', lw=0.5)
ax[1].axvline(f2, c='k', lw=0.5)
ax[1].axhline(phi1, c='k', lw=0.5)
ax[1].axhline(phi2, c='k', lw=0.5)
ax[1].legend(loc=(1.01, 0));
../_images/370edad5058c31ee4de103f21ecefb467a1298a4039d598032d2e2080a5ec3ad.png

Primerjava izbranih faznih zamikov signala in vrednosti faznega spektra:

print(phi1, phi2)
print(np.angle(S[i1]), np.angle(S[i2]))

1.5707963267948966 0.7853981633974483
1.5707963267948966 0.7853981633974466
Contents