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 (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. Demonstrirajte 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/e67c017ff791755d8caa72dccc8ed3741bf95d6a8e53311eeb199bbaee5769b8.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/e1c7c544621db9c432dee87d11b027fe46e756d733f95ece5bfa4490cd6ea575.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/998449e3dba038465ea28f500d1694b3e8cc2ef69bc9fadd2df27759eab584d3.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/efcd54ec541ebb9a75d1458cb72d7c059a1d4e5864156e25484d7c3fd0ddb40a.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/e6acb190558e6b4dbec13278ab93edf7784ab9970fcd90ec7aef5b48e30b772b.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/4827386ae63b89f732f7af71e28f9c569759a94c2b602cfcd232905c0373afee.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/137299fdf476c98f0e920259a223bd89c0d2c6a450c7199f69792cdd78f3a890.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/fbf640e4bab5b8ea23d37ff3b412d3d2e41873cb2ec114dc69a8414c5fc14fae.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/65ea90ad35c58889323776d1f955e1425bfe50b7ec451aebecd92477f9dfd945.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/719d6e0b7a9ee9080349b60939b299891c8f404a8dac2aafdda7cbc535323aa3.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/3e2a359557dbfcb1127c1874296feb6fcd9c0a726567346feb22a3d15ad3b75c.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/29eafe3f3cf4cd50e4700722dc53da6f4eaf02d9cd215c1682319bf39aa78e82.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/3896cb4a02a7a3259b3c1995e9b6a8d019c87b671761f9d846eac80e8854fc9d.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 0x7f4e27773110>
../_images/b5a8f7a3335c5ecb197e7f45716c48ac3f60bac33076fa45fac822b3e8cb06d8.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/005a96ec8ac68f2b327f62ff88d31596f875c583c12120d4831e1434737f7e61.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/07ee080b2127d93e3199baab2d1b0747190d9a6884826d158323e725be313cfc.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/2750fb665d8408bd44c25e1fe8bde10c3520ab7401dbaec360f765699cff1088.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/155fd0636f59ff061caa96658e793b3e6022fa1eee976117619024401ba6432e.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/aa1829952c461fa239104b41ed65bc0667b119d2874d781deb0a1e4ba311dcae.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/7a2378ebc0e04775cbdd738c3d8f53615fa9904d4d28cd87d0c0b408fafc4665.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/70a4e1b8cf80091886b9abe2262a4e83d8f7952528a1402d0a0a92234552de38.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/c983fa659b6b30798c5b5b670da396f784be8a5c0ad25db7847df0cb137fef59.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/a59ebd57f589c53b8635b392499a856cab320a7d8cc8d6cb156f4e429d01db32.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/4c6401a86de0e63263b1f6b27b9c9fa26159eeaeff9303720c83bbd87fd6170b.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/a5d6bd6abf4fe1b9a04b4744e695e3a456a43abee28507bdc3b72abc6098b847.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/c47d85557da5654f145fef4b9f489a9b6d35383a384b11a0a962bd3a4998d1ad.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/038601837b76bbc298590ba5b98da115b7a0f59867fd3b4f23c78e5dcd64f2ac.png

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

Contents