Ctrl+K

Lastnosti Fourierove transformacije

Hide code cell source 
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

Lastnosti Fourierove transformacije#

Na tej vaji si bomo iz praktičnega vidika pogledali lastnosti Fourierove transformacije, s katerimi imamo pogosto opravka pri merjenju signalov:

  • linearnost,

  • časovno skaliranje,

  • (časovni obrat),

  • časovni premik,

  • (amplitudna modulacija - frekvenčni premik),

  • odvajanje in integriranje,

  • konvolucija in produkt funkcij.

Naloga#

Domača naloga

Z uporabo generatorja signalov in zajemnega sistema Arduino pripravite signale, s katerimi boste demonstrirali lastnosti Fourierove transformacije, določene v podatkih naloge.

Priprave in zajema signalov se lotite v skupini s kolegi, ki imajo enake podatke. LabView program za zajem lahko prenesete v obliki zip arhiva.

Pripravite kratko poročilo v okolju Jupyter Notebook (od 3 do 10 celic s kodo), v katerem naj bodo razvidni podatki naloge, ter demonstriran vpliv izbrane lastnosti Fourierove transformacije na zajeta signala (zglede poiščite v predlogah za predavanja in vaje).

Poročilo oddajte v .pdf obliki (glejte navodila za oddajo domačih nalog).

Dodatek: Raziščite lastnost konvolucije / množenja signalov v povezavi s Fourierovo transformacijo in jo prikažite ter komentirajte na primeru uporabe oken (uporabite scipy.signal.windows) z enim izmed signalov, zajetih na vaji.

Linearnost#

Linearnost

\[ \mathcal{F}\{a\,x(t)+b\,y(t)\} = a\,X(f) + b\,Y(f). \]
fs = 100 # Hz
dt = 1 / fs
t = np.arange(0, 2, dt)
A = 5
p = 4
o = 2

x = A * np.sin(2*np.pi*t * p) + o*2
y = 2*A * signal.square(2*np.pi*t * p/2) + o
Hide code cell source 
plt.figure()
plt.plot(t, x, label='sin')
plt.plot(t, y, label='square')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/2fe1fcfaab5ce9e3ce3863b838421909a5799482b81f30524b996d1b8924b62e.png 
X = np.fft.rfft(x) / len(t) * 2
Y = np.fft.rfft(y) / len(t) * 2
freq = np.fft.rfftfreq(len(t), dt)
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='X')
plt.plot(freq, np.abs(Y), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.legend(loc=(1.01, 0));
../_images/e675afb08d82f11ea87eb717c48fb81b11bb2cc3d04e8a4b9645702e035431ed.png

Prikaz lastnosti linearnosti:

v = x + y
V = np.fft.rfft(v) / len(t) * 2
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(V), label='$\\mathcal{F}(x+y)$')
plt.plot(freq, np.abs(X + Y), label='$\\mathcal{F}(x) + \\mathcal{F}(y)$')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.legend(loc=(1.01, 0));
../_images/d3516d4a64f8717b5ecb3b7b2f4fb4cc5624a58332ee2ad184f4b4845593229f.png

Uporaba linearnosti pri korekciji naloženega signala:

V_Y = V - Y
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(V_Y), label='$\\mathcal{F}(x+y) - \\mathcal{F}(y)$')
plt.plot(freq, np.abs(X), label='$\\mathcal{F}(x)$')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.legend(loc=(1.01, 0));
../_images/65a80d6957688c0f79c1863d8b280ccc7d65e6c9e4a6305a0b0ad02c45479bc4.png

Inverzna FFT razlike vsote signalov in naloženega signala:

v_y = np.fft.irfft(V_Y/2*len(t))
Hide code cell source 
plt.figure()
plt.plot(t, x, label='x')
plt.plot(t, v_y, label='$\\mathrm{IFFT}(\\mathcal{F}(x+y) - \\mathcal{F}(y))$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/9906e94d504d6013ef4eaa20c7bf6a7c56981fc45f86fa76d77b4b0082780a97.png

Časovno skaliranje#

Časovno skaliranje

\[ \mathcal{F}\{x(a\,t)\} = \frac{1}{|a|}\,X(f/a). \]
A = 1
p = 2
o = 0
a = 2

sin = lambda t, A, p, o: A*np.sin(2*np.pi*t*p) + o
x = sin(t, A, p, o)

at = np.arange(0, 2, a*dt)
y = sin(at, A, p, o)
Hide code cell source 
plt.figure()
plt.plot(t, x, label='x(t)')
plt.plot(at, y, label='x(at)')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/a104d7a90c13027228ef84d638944bed885d021823effc4179b5827c8714fae3.png

Trajanje meritve je enako, a pri časovno skaliranem koraku \(a \, \Delta t\):

Hide code cell source 
plt.figure()
plt.plot(x, label='x(t)')
plt.plot(y, label='x(at)')
plt.xlabel('vzorec [/]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/91afaa02619e9a17d28fbb9af2a7c223edac39e74681f37334fb23f405ca51d4.png 
X = np.fft.rfft(x)
Y = np.fft.rfft(y)
freq = np.fft.rfftfreq(len(t), dt)
freq_a = np.fft.rfftfreq(len(at), a*dt)

Če časovno skaliranje poznamo in upoštevamo pri pretvorbi je rezultat pravilen:

Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X/len(t)*2), label='X')
plt.plot(freq_a, np.abs(Y/len(at)*2), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.legend(loc=(1.01, 0));
../_images/533984b0b57d59867544b08ccd64cf70d775e04bffe3a346bd83a35f56117b5a.png

Časovno skaliranje nam v tem primeru povzroča težave, če se ga pri meritvi ne zavedamo:

freq_a_s = np.fft.rfftfreq(len(at), dt)
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='X')
plt.plot(freq_a_s, np.abs(Y), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/3d0e393e270dd7bc0c2b0784bb592fbc4d950a3e933df130326d9a3b964fadd0.png

Korekcija časovnega skaliranja

\[ \mathcal{F}\{x(a\,t)\} = \frac{1}{|a|}\,X(f/a). \]

Signal smo pretvorili v frekvenčno domeno pri frekvencah \(f/a\), sedaj ga še narišimo pri \(f/a\):

freq_a_popravek = freq_a_s / a

Popravek amplitude:

Y_popravek = np.abs(a) * Y
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='X')
plt.plot(freq_a_popravek, np.abs(Y_popravek), label='Y_popravljen')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/bac8c34a2846c3f21db9a6bdd23ccad2d36ec7b2605f022a27701720f973ce85.png

Časovni premik#

Časovni premik

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

square = lambda t, A, p, o: A*signal.square(2*np.pi*t*p) + o

x = square(t, A, p, o)
y = square(t-t_0, A, p, o)
Hide code cell source 
plt.figure()
plt.plot(t, x, label='x(t)')
plt.plot(t, y, label='x(t - t_0)')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/b420fe6518842e9de947857f1f9d75da5a4dd87e60acde52f032d8fa3937b256.png 
X = np.fft.rfft(x)
Y = np.fft.rfft(y)
freq = np.fft.rfftfreq(len(t), dt)
Hide code cell source 
fig, ax = plt.subplots(2, 1)
ax[0].plot(freq, np.abs(X/len(t)*2), label='X')
ax[0].plot(freq, np.abs(Y/len(t)*2), label='Y')
ax[0].set_ylabel('$|\\mathcal{F}(x)|$ [/]')
ax[0].set_xlim(0, 20)
ax[0].legend(loc=(1.01, 0));

ax[1].plot(freq, np.angle(X/len(t)*2), label='X')
ax[1].plot(freq, np.angle(Y/len(t)*2), label='Y')
ax[1].set_xlabel('f [Hz]')
ax[1].set_ylabel('$\\angle \\mathcal{F}(x)$ [rad]')
ax[1].set_xlim(0, 20)
ax[1].legend(loc=(1.01, 0));
../_images/9edf0072e76ade0a42718a79c902406ad6acaeca1aa4c7a331cf2e7224d1e2cc.png 
X_shift = np.exp(-1j*2*np.pi*freq*t_0) * X
Hide code cell source 
fig, ax = plt.subplots(2, 1)
ax[0].plot(freq, np.abs(X/len(t)*2), label='X')
ax[0].plot(freq, np.abs(X_shift/len(t)*2), label='shifted X')
ax[0].set_ylabel('$|\\mathcal{F}(x)|$ [/]')
ax[0].set_xlim(0, 20)
ax[0].legend(loc=(1.01, 0));

ax[1].plot(freq, np.angle(X/len(t)*2), label='X')
ax[1].plot(freq, np.angle(X_shift/len(t)*2), label='shifted X')
ax[1].set_xlabel('f [Hz]')
ax[1].set_ylabel('$\\angle \\mathcal{F}(x)$ [rad]')
ax[1].set_xlim(0, 20)
ax[1].legend(loc=(1.01, 0));
../_images/d32f99a9e1b839e45e0f8aa38c21cd2bb00e0814cb91e92fe7914f6e97bbce61.png

Uspešnost časovnega premika preverimo v časovni domeni:

x_shift = np.fft.irfft(X_shift)
Hide code cell source 
plt.figure()
plt.plot(t, x_shift, label='$\\mathrm{IFFT}(X_{\\mathrm{shift}})$')
plt.plot(t, y, label='x(t - t_0)')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0))

<matplotlib.legend.Legend at 0x18ba7386fc0>
../_images/42b64408269e730f6bef434fb907e1863611ba5fc2cb2039422b59cc8e01b81b.png

Amplitudna modulacija#

Amplitudna modulacija (frekvenčni zamik)

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

oziroma:

\[ \mathcal{F}\{x(t)\,\cos(2\pi\,f_0\,t)\} = \frac{1}{2}\big(X(f-f_0)+X(f+f_0)\big). \]
A = 1
p = 2
o = 0
f_0 = 1.

x = sin(t, A, p, o)
mod = np.cos(2*np.pi*f_0*t)
y = x * mod
Hide code cell source 
plt.figure()
plt.plot(t, x, label='$x(t)$')
plt.plot(t, mod, label='modulacijski signal')
plt.plot(t, y, label='$x_{\\mathrm{moduliran}}(t)$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/166714df0fa825c7d53dee63c982fa9dafc394aead5bc37e8f6778dddf60efbb.png 
X = np.fft.rfft(x)
Y = np.fft.rfft(y)
freq = np.fft.rfftfreq(len(t), dt)
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X/len(t)*2), label='X')
plt.plot(freq, np.abs(Y/len(t)*2), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/3cc9d7403a1a2606da89aa8425ba39b0480b6e62bb9d6ca9ae05b5ed33f11d93.png

Odvajanje in integriranje#

Odvajanje in integriranje

\[ \mathcal{F}\{\dot x(t)\} = \textrm{i}\,2\pi\,f\,X(f) \]

in obratno:

\[ \mathcal{F}\big\{\int_{-\infty}^t x(\tau)\,\textrm{d}\tau\big\} = \frac{X(f)}{\textrm{i}\,2\pi\,f}. \]
A = 1
p = 2
o = 1.5

x = sin(t, A, p, o)

Spomnimo se:

\(\frac{\mathrm{d}}{\mathrm{d} \, t}A \, \sin(a \, x) + o = A \, a \, \cos(a \, x)\)

y = A*2*np.pi*p*np.cos(2*np.pi*t*p)

np.allclose(y, np.gradient(x, dt, edge_order=2), rtol=1e-2)

True
Hide code cell source 
plt.figure()
plt.plot(t, x, label='$x(t)$')
plt.plot(t, y, label='$y = \\dot{x}(t)$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/49009a6f3350577d53947c329fe735f8aaaaea9268c377902538633ccfca3bd9.png 
X = np.fft.rfft(x) / len(t)
Y = np.fft.rfft(y) / len(t)
X[1:] *= 2
Y[1:] *= 2
freq = np.fft.rfftfreq(len(t), dt)
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='X')
plt.plot(freq, np.abs(Y), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/5b0609a2c98daab5c295e080e7fbc38c0176b99b8918dbacc406532d352521b4.png

Odvod v frekvenčni domeni:

\(\mathcal{F}(\dot{x}(t)) = \mathrm{i} \, \omega \, \mathcal{F}(x(t))\)

omega = 2*np.pi*freq
dX = 1j*omega * X
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(dX), label='$\\mathrm{i} ~ \\omega X$')
plt.plot(freq, np.abs(Y), label='Y')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/d8bd2985701a59abfd96d9e6f68104bd7e550fb645461a839d9a0769ea56a17d.png

Integral v frekvenčni domeni:

\(\mathcal{F}(x(t)) = 1 / (\mathrm{i} \, \omega) \cdot \mathcal{F}(\dot{x}(t))\)

omega[0] = 1e-6
ing_Y = 1 / (1j*omega) * Y
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='$X$')
plt.plot(freq, np.abs(ing_Y), label='$1/(\\mathrm{i} ~ \\omega) ~ Y$')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 10)
plt.legend(loc=(1.01, 0));
../_images/82334cc953ad4232eadb0f7d7759163999b09a5d5faee341ed12ab00f0add3e3.png

Konvolucija in produkt funkcij#

Konvolucija in produkt funkcij

\[ \mathcal{F} \big\{ x(t)*y(t) \big\} = X(f) \, Y(f) \]

in

\[ \mathcal{F} \big\{ x(t) \, y(t) \big\} = X(f)*Y(f) \]

Note

Znak \(*\) označuje konvolucijo:

\[ x(t)*y(t) = \int_{-\infty}^{+\infty} x(\tau)\,y(t-\tau)\,\textrm{d}\tau \]

Ker je produkt v frekvenčni domeni numerično veliko učinkovitejša operacija od konvolucije v časovni domeni, je to ena najpomembnejših lastnosti Fourierove transformacije!

Vizualizacija konvolucije na diskretnem signalu:

convolve

A = 1
p = 3
o = 1.5

x = square(t, A, p, o)

X = np.fft.rfft(x) / len(t)
X[1:] *= 2
freq = np.fft.rfftfreq(len(t), dt)
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='X')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 20)
plt.legend(loc=(1.01, 0));
../_images/e6396962974a3f1f32d4e58101d5a44bd97237a5d515bdc9051cb7c91b516c8c.png

Pripravimo pasovni frekvenčni filter, ki prepušča le v območju \((1.5, 4.5)\) Hz:

freq_band = np.array([1.5, 4.5]) 
b, a = signal.butter(3, freq_band, btype='bandpass', fs=fs)

Pripravimo frekvenčno prenosno funkcijo filtra \(H(f)\):

w, H = signal.freqz(b, a, len(freq), whole=False, fs=fs)
Hide code cell source 
plt.figure()
plt.plot(w, np.abs(H), label='filter')
plt.plot(freq, np.abs(X), label='X')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 20)
plt.legend(loc=(1.01, 0));
../_images/3536bafac3fa92c2c34c7f85a5ab01d97965ca7397f815e1e71003b4bf322ea8.png

Filtrirajmo signal z množenjem v frekvenčni domeni:

X_filtriran = X * H
Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X), label='$X$')
plt.plot(freq, np.abs(X_filtriran), label='$X \\cdot H$')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.xlim(0, 20)
plt.legend(loc=(1.01, 0));
../_images/f36ac887a54b1cea3c267ab0777e664a0f4c887e2ba4b46be42365e5dae00c1b.png

Poglejmo še impulzno prenosno funkcijo filtra \(h(t)\) (odziv filtra v času):

h = np.fft.irfft(H)
Hide code cell source 
plt.figure()
plt.plot(t, x, label='$x(t)$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=2);
ax2 = plt.gca().twinx()
ax2.plot(t, h, 'C1', label='filter')
ax2.set_ylabel('filter [/]')
ax2.legend(loc=1);
../_images/919e3b7fe1fb5f60a804831cfc723fa527c502a3e9954d60257015e757153d8d.png

Filtrirajmo signal s konvolucijo v časovni domeni:

x_filtriran = signal.convolve(x, h, mode='full')[:len(x)]
Hide code cell source 
plt.figure()
plt.plot(t, x, label='$x(t)$')
plt.plot(t, x_filtriran, label='$x * h$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.legend(loc=(1.01, 0));
../_images/dff05d8660ba978ed3fe93d409fefac7cf541cee9822680e53aa51d55dc13073.png

Primerjajmo rezultata:

Hide code cell source 
plt.figure()
plt.plot(freq, np.abs(X_filtriran), label='$X \\cdot H$')
plt.plot(freq, np.abs(np.fft.rfft(x_filtriran)*2/len(t)), label='$\\mathcal{F}(x*h)$')
plt.xlabel('f [Hz]')
plt.ylabel('$|\\mathcal{F}(x)|$ [/]')
plt.title('Frekvenčna domena')
plt.xlim(0, 20)
plt.legend(loc=(1.01, 0));
../_images/b9e62e65825caafe133bb365fc5a9287cb2e18a613e1b8a888cf97093cac55be.png
Hide code cell source 
plt.figure()
plt.plot(t, x_filtriran, label='$x * h$')
plt.plot(t, np.fft.irfft(X_filtriran*len(t)/2), label='$\\mathrm{IFFT}(X \\cdot h)$')
plt.xlabel('t [s]')
plt.ylabel('x [/]')
plt.title('Časovna domena')
plt.legend(loc=(1.01, 0));
../_images/572a5248ef3e75fe78bc408cbcaa2050378835f959c0eb0a5ee51c8a659c4f35.png

Razlika je posledica implementacije funkcije konvolucije (scipy.signal.convolve). Preverite vpliv parametra mode!

Contents