Show code cell source
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = [12.0, 7.0]
mpl.rcParams['font.size'] = 14
import warnings
warnings.filterwarnings('ignore')
Linearni, časovno invariantni (LTI) sistemi 2 - prenosne funkcije#
Domača naloga#
Domača naloga
N podlagi meritev nosilca, ki ste jih v skupini zajeli na prejšnji vaji, izračunajte in prikažite frekvenčne prenosne funkcije pospešenosti za vse pomerjene lokacije vzbujanja na nosilcu. Pri obdelavi zajetih signalov ustrezno uporabite eksponentno okno ter “force” okno s parametri iz podatkov naloge.
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: Na podlagi pripravljenih FPF pri različnih točkah nosilca ocenite prvo lastno obliko obravnavanega nosilca.
Lastna oblika podaja razmerje amplitud odziva sistema na vzbujanje v različnih točkah \(i\):
kjer so \(A^r_i\) t. i. modalne konstante sistema pri lastni frekvenci \(\omega_r\).
Modalne konstante lahko ocenite s poljubno metodo (glejte spodnji zaznamek).
Za identifikacijo modalnih konstant lahko uporabite tudi paket sdypy-EMA. Pomagate si lahko s temle primerom uporabe.
Metoda največjih amplitud
Ena najpreprostejših metod modalne identifikacije je metoda največjih amplitud (Peak amplitude, tudi Peak Picking method):
\(\omega_r = \text{arg max}(|\alpha(\omega)|)\) (krožna frekvenca resonančnega vrha)
\(\eta_r = \frac{\omega_a - \omega_b}{\omega_r}\) (razmernik dušenja), kjer sta \(\omega_a\) in \(\omega_b\) t. i. točki polovične moči v okolici \(\omega_r\):
\(|\alpha(\omega_{a, b})| = \frac{|\alpha(\omega_r)|} {\sqrt{2}}\)
Iz FPF za posamezno obravnavano točko \(i\) pa lahko določite modalne konstante:
\(A^r_i = |\alpha_i(\omega_r)| \, \omega_r^2 \, \eta_r\)
Show code cell source
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import lvm_read
Frekvenčna prenosna funkcija \(\alpha(\omega)\) predstavlja razmerje med (kompleksnimi) amplitudami harmonskega odziva in vzbujanja opazovanega sistema:
Impulzna prenosna funkcija \(h(t)\) podaja odziv sistema na impulzno motnjo.
Odziv na poljubno vzbujanje (signal \(f(t)\)) določa naslednja zveza:
kjer \(*\) označuje konvolucijo.
Note
Impulzna prenosna funkcija in frekvenčna prenosna funkcija sta Fourierov par:
Opis LTI sistemov#
Odzivni model LTI sistema opišemo z njegovo prenosno funkcijo.
Poglejmo primer enostavnega linearnega oscilatorja, ki ga opisuje gibalna enačba:
Frekvenčna domena#
Frekvenčna prenosna funkcija dušenega linearnega oscilatorja je podana z:
Laplace-ova domena
Fourierova transformacija je posebna oblika Laplace-ove transfomacije, ki jo pogosto uporabimo za pretvorbo diferencialnih enačb v algebrajske enačbe.
Velja naslednja zveza med frekvenčno (\(\omega\)) in Laplace-ovo (\(s\)) domeno:
\(s = \mathrm{i} \, \omega\)
Zgornjo prenosno funkcijo linearnega oscilatorja v \(s\)-domeni tako zapišemo:
Poglejmo konkreten primer:
m = 2
k = 50
c = 2.5
freq = np.linspace(0, 10, 101)
w = 2*np.pi*freq
def H_f(freq, m, k, c):
"""
Frekvenčna prenosna funkcija linearnega oscilatorja.
"""
omega = 2*np.pi*freq
return 1 / (-omega**2*m + 1j*omega*c + k)
def H_s(freq, m, k, c):
"""
Prenosna funkcija linearnega oscilatorja v s-domeni.
"""
omega = 2*np.pi*freq
s = 1j*omega
return 1 / (m*s**2 + c*s + k)
Show code cell source
plt.plot(freq, np.abs(H_f(freq, m, k, c)), label=r'$\alpha({\omega})$')
plt.plot(freq, np.abs(H_s(freq, m, k, c)), '--', label='$H(s)$')
plt.axvline(x=np.sqrt(k/m) / (2*np.pi), c='k', lw=0.5)
plt.legend()
plt.xlabel('f [Hz]')
plt.ylabel('|H(f)|');

Naloga 1 (10 minut): uporaba scipy.signal.lti
Sistem linearnega oscilatorja definirajte z uporabo scipy.signal.lti
. Pomagajte si z dokumentacijo.
Pri tem uporabite prenosno funkcijo sistema v \(s\)-domeni:
Primerjajte frekvenčno in impulzno prenosno funkcijo, pripravljeni na podlagi zgornje enačbe, s tistima, ki ju dobite z lti.freqresp
ter lti.impulse
.
T = 1 / (freq[1] - freq[0])
dt = 1/ (2*freq[-1])
t = np.arange(0, T+dt, dt)
Spomnimo se, frekvenčna prenosna funkcija \(\alpha(\omega)\) in impulzna prenosna funkcija \(h(t)\) sta Fourierov par:
h = np.fft.irfft(H_f(freq, m, k, c), n=2*(len(freq)-1)+1)
Show code cell source
plt.plot(t, h)
plt.xlabel('t [s]')
plt.ylabel('$h(t)$');


Opazimo, da sta izračunani impulzni prenosni funkciji različno normirani.

Naloga 2 (20 minut): lastnosti LTI sistemov
Z uporabo zgoraj pripravljene impulzne prenosne funkcije \(h(t)\) ter lti.output
na dva načina prikažite načelo superpozicije:
in lastnost časovne invariantnosti:
LTI sistema.
Uporabite spodaj pripravljena signala \(f_1\), \(f_2\) ter vrednosti konstant \(a_1\), \(a_2\).
a_1 = 25
a_2 = -35
f_1 = np.zeros_like(t)
f_1[10:20] = np.sin(2*np.pi*np.arange(10)/20)
f_2 = np.zeros_like(t)
f_2[32:36] = np.sin(2*np.pi*np.arange(4)/8)
Show code cell source
plt.plot(t, a_1*f_1, label='$a_1 \, f_1(t)$')
plt.plot(t, a_2*f_2, label='$a_2 \, f_2(t)$')
plt.xlabel('t [s]')
plt.ylabel('$f(t)$')
plt.legend();



Obdelava signalov za oceno impulzne prenosne funkcije#
Naloga 3 (30 minut): Prenosne funkcije nosilca
Na primeru nosilca, obravnavanega na prejšnji vaji, ocenimo izmerjene frekvenčne prenosne funkcije nosilca.
Naložimo datoteko z meritvami:
meritev_pot = 'data/06/meritev.lvm'
meritev = lvm_read.read(meritev_pot)
Izberimo eno od obravnavanih točk.
Za demonstracijo pomena obdelave zajetih signalov uporabimo le del zajetih časovnih vrst (\(T_{\text{analiza}}\)).
točka = 7
T_analiza = 0.5
dt = meritev[točka]['Delta_X'][0]
N = int(T_analiza * 1/dt)
print(meritev[točka]['Y_Unit_Label'])
['Newtons', 'g', '']
f, a_g = meritev[točka]['data'][:N].T
a = a_g * 9.81
t = np.arange(len(f)) * dt
Show code cell source
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(t, a, alpha=1, label='$a(t)$')
ax2.plot(t, f, alpha=1, c='C1', label='$f(t)$')
ax.legend(loc=1)
ax2.legend(loc=2)
ax2.set_ylabel('F [N]')
ax.set_ylabel('a [m/s$^2$]')
ax.set_xlabel('t [s]');

Opazimo, da se v opazovanem časovnem intervalu zajet odziv sistema ne izniha popolnoma.
Poglejmo amplitudna spektra zajetih signalov za izbrano točko:
F = np.fft.rfft(f) / len(f)
F[1:] *= 2
A = np.fft.rfft(a) / len(a)
A[1:] *= 2
freq = np.fft.rfftfreq(len(t), dt)
Show code cell source
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.semilogy(freq, np.abs(A), alpha=1, label='$A(f)$')
ax2.semilogy(freq, np.abs(F), alpha=1, c='C1', label='$F(f)$')
ax.legend(loc=1)
ax2.legend(loc=2)
ax2.set_ylabel('|F [N]|')
ax.set_ylabel('a [m/s$^2$]')
ax.set_xlabel('f [Hz]');

Ocenimo frekvenčno prenosno funkcijo:
FRF = A / F
Show code cell source
plt.semilogy(freq, np.abs(FRF))
plt.ylabel('$|H(f)|$')
plt.xlabel('f [Hz]');

Ekponentno okno#
Da zadostimo predpostavki periodičnosti signala v primeru, ko se odziv sistema v zajetem časovnem intervalu ne izniha, uporabimo t. i. eksponentno okno:
Namesto parametra časovne konstante \(\tau\) pogosto želimo podati končno amplitudo \(p\) (v obliki deleža začetne amplitude), ki jo okno doseže pri času \(T\). Velja:
p = 0.01 # 1 %
tau = -(len(t) - 1) / np.log(p)
w_exponential = signal.windows.exponential(len(t), center=0, tau=tau, sym=False)
w_exponential[-1]
0.010000000000000004
Primer uporabe eksponentnega okna na numeričnem modelu#
k2 = (100*2*np.pi)**2
m2 = 1
c2 = 5.
sistem_primer = signal.lti([1], [m2, c2, k2])
f_primer = np.zeros_like(t)
f_primer[1000:1010] = 1.
t, x_primer, _ = sistem_primer.output(f_primer, t)
Show code cell source
fig, ax = plt.subplots()
ax.plot(t, x_primer, label='odziv')
ax2 = ax.twinx()
ax2.plot(t, f_primer, 'C1', label='vzbujanje');
ax.legend(loc=1)
ax2.legend(loc=2)
ax2.set_ylabel('F [N]')
ax.set_ylabel('a [m/s$^2$]')
ax.set_xlabel('t [s]');

X_primer = np.fft.rfft(x_primer) / len(t)
X_primer[1:] *= 2
Show code cell source
plt.semilogy(freq, np.abs(X_primer))
plt.xlim(0, 500);
plt.ylabel('$|H(f)|$')
plt.xlabel('f [Hz]');

Ker v opazovanem časovnem oknu nismo zajeli popolnoma periodičnega signala, v amplitudnem spektru opazimo pojav frekvenčnega odtekanja (ang. Spectral Leakage).
Odpravimo ga lahko z ustrezno uporabo oken.
X_okno = np.fft.rfft(x_primer*w_exponential) / len(t)
X_okno[1:] *= 2
Show code cell source
plt.semilogy(freq, np.abs(X_primer), label='brez okna')
plt.semilogy(freq, np.abs(X_okno), label='z eksponentnim oknom')
plt.xlim(0, 500);
plt.legend()
plt.ylabel('$|H(f)|$')
plt.xlabel('f [Hz]');

F_primer = np.fft.rfft(f_primer) / len(t)
F_primer[1:] *= 2
w, H_primer = sistem_primer.freqresp(2*np.pi*freq)
H_okno = X_okno/F_primer
Show code cell source
plt.semilogy(freq, np.abs(H_okno), label='z eksponentnim oknom')
plt.semilogy(freq, np.abs(H_primer), label='sympy.lti')
plt.xlim(0, 500);
plt.ylim(1e-9, 1e-2)
plt.legend()
plt.ylabel('$|H(f)|$')
plt.xlabel('f [Hz]');

Uporaba eksponentnega okna na meritvi nosilca#
Show code cell source
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(t, a, alpha=1, label='$a(t)$')
ax2.plot(t, w_exponential, alpha=1, c='C1', label='$w(t)$')
ax2.set_ylim(-1.1, 1.1)
ax.legend(loc=1)
ax2.legend(loc=2)
ax2.set_ylabel('w(t) [/]')
ax.set_ylabel('a [m/s$^2$]')
ax.set_xlabel('t [s]');

a_w = a * w_exponential
A_w = np.fft.rfft(a_w) / len(a_w)
A_w[1:] *= 2
Show code cell source
fig, ax = plt.subplots()
ax.semilogy(freq, np.abs(A), alpha=1, c='C1', label='$\mathcal{F}(a)$')
ax.semilogy(freq, np.abs(A_w), alpha=1, label='$\mathcal{F}(a \cdot w)$')
ax.legend(loc=1)
ax.set_ylabel('A(f) [N]')
ax.set_xlabel('f [Hz]');

Force okno#
Da zmanjšamo vpliv šuma na vhodnem signalu \(f(t)\) uporabimo t. i. “force” okno, s katerim izničimo vrednosti signala sile po koncu udarca.
Parameter “force” okna je dolžina segmenta okna (v obliki deleža skupne dolžine signala) z enotsko amplitudo, ki jo bomo označili z \(l_w\).
Pri pripravi “force” okna si lahko pomagamo z oknom kosinusnega upada scipy.signal.windows.tukey
.
l_w = 0.1
N_w = int(l_w * len(t))
w_force = np.zeros_like(t)
w_force[:N_w] = signal.windows.tukey(len(t)*2, alpha=0.05)[-N_w:]
Ker smo na signalu odziva uporabili eksponentno okno, moramo eksponentno okno uporabiti tudi na signalu sile.
w_f = w_force * w_exponential
Show code cell source
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(t, w_force, alpha=1, c='C2', label='$w(t)$')
ax.plot(t, w_f, alpha=1, c='C4', label='$w_f(t) \cdot w_e(t)$')
ax2.plot(t, f, alpha=1, c='C1', label='$f(t)$')
ax.legend(loc=0)
ax.set_ylabel('w(t) [/]')
ax2.set_ylabel('f(t) [N]')
ax.set_xlabel('t [s]');

Uporaba force okna:
f_w = f * w_f
F_w = np.fft.rfft(f_w) / len(f_w)
F_w[1:] *= 2
Show code cell source
fig, ax = plt.subplots()
ax.semilogy(freq, np.abs(F), alpha=1, c='C1', label='$\mathcal{F}(f)$')
ax.semilogy(freq, np.abs(F_w), alpha=1, lw=2, label='$\mathcal{F}(f \cdot w)$')
ax.legend(loc=1)
ax.set_ylabel('F [N]')
ax.set_xlabel('f [Hz]');

FRF z okni#
FRF_w = A_w / F_w
Show code cell source
plt.semilogy(freq, np.abs(FRF), label='brez oken')
plt.semilogy(freq, np.abs(FRF_w), label='z okni')
plt.ylabel('$|H(f)|$')
plt.xlabel('f [Hz]')
plt.legend()
plt.xlim(0, 5000);
