Show code cell source
import pickle
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
Naključno vzbujanje - ocena frekvenčnih prenosnih funkcij#
Spomnimo se:
Note
Frekvenčna prenosna funkcija
Določamo razmerje med odzivom sistema
Namen te vaje je spoznati metode procesiranja signalov za določitev cenilk frekvenčnih prenosnih funkcij linearnega sistema na podlagi tovrstnih meritev.
Domača naloga#
Domača naloga
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. Vsebuje naj rezultate in kratek komentar izvedenih meritev ter izračunanih cenilk frekvenčne prenosne funkcije glede na podatke naloge. Poročilo naj vsebuje:
prikaz PSD vaše meritve (povprečje vseh segmentov, v točki, določeni v podatkih naloge 9),
prikaz frekvenčne prenosne funkcije v vaši izmerjeni točki (
), ki jo ocenite z .prikaz frekvenčne prenosne funkcije v vaši izmerjeni točki (
), izračunane s cenilko ali (glede na podatke naloge). Pri tem povprečite spektre ( ) izmerjenih signalov, razdeljenih na segmente dolžine 2 s.
Dodatek:
Izračunajte in prikažite obe cenilki frekvenčne prenosne funkcije (
, ).Križne- in avtospektre signalov pri tem izračunajte s povprečenjem prekritih segmentov in uporabo oknjenja,
Dolžina posameznega segmenta izračuna PSD in CSD je določena v podatkih naloge, uporabite 50% prekrivanjem segmentov,
Na posameznem segmentu za zmanjševanje pojava frekvenčnega odtekanja uporabite okno, določeno v podatkih naloge.
Izračunajte in prikažite tudi koherenco meritve FRF,
.
Poročilo oddajte v .pdf
obliki (glejte navodila za oddajo domačih nalog).
Branje izmerjenih signalov#
Parametri izvedene meritve
Frekvenčne prenosne funkcije jeklenega nosilca bomo izračunali na podlagi meritev, izvedenih pri vaji 9.
Pri procesiranju zajetih signalov so pomembni tudi parametri izvedene meritve. Če ste meritvah sledili navodilom, so osnovni parametri vaših signalov naslednji (če ste pri meritvi uporabili drugačne nastavitve, vrednosti po potrebi prilagodite!):
Osnovni parametri zajema:
Dva kanala: sila in pospešek (IEPE).
Frekvenca vzorčenja: 25600 Hz,
Trajanje vsakega zajetega segmenta: 2 s,
Število zajetih segmentov: 5.
Osnovne značilnosti meritve:
11 označenih merilnih mest (na pozicijah
mm, razmik mm).Vzbujanje na lokaciji
, meritve odziva na lokacijah .Trajanje vsakega generiranega segmenta: 2 s,
Profil naključnega vzbujevalnega signala konstantne amplitude na območju
,Povprečne amplitude vzbujevalne sile: 1 N RMS.
Branje shranjene datoteke - meritev v eni izmed označenih točk
pot = './data/10/meritev.pkl'
meritev = pickle.load(open(pot, 'rb'))
print(meritev.keys())
dict_keys(['data', 'fs', 'channels'])
Frekvenca vzorčenja in shranjeni kanali
fs = meritev['fs']
print(fs, meritev['channels'])
25600.0 ['Force', 'Acceleration']
Oblika shranjenih signalov: (segmenti, kanali, čas)
print(meritev['data'].shape)
(5, 2, 51200)
Izbira segmenta za analizo
force, acc = meritev['data'].transpose(1, 0, 2) # transponiranje signalov tako, da imamo na prvi osi oba pomerjena kanala, silo in pospešek
force.shape
(5, 51200)
time = np.arange(force.shape[1]) / fs
F = np.fft.rfft(force, axis=-1) / force.shape[-1] * 2 # izračun FFT: po časovni osi
A = np.fft.rfft(acc, axis=-1) / acc.shape[-1] * 2
F.shape
(5, 25601)
freq = np.fft.rfftfreq(len(time), 1/fs)
freq.shape
(25601,)
Izris amplitudnih spektrov
Show code cell source
segment = 0
fig, ax = plt.subplots(2, 1)
plt.suptitle(f'Segment [{segment}]')
ax[0].semilogy(freq, np.abs(A[segment]), label='pospešek')
ax[0].semilogy(freq, np.abs(F[segment]), label='sila')
ax[0].set_xlabel('frekvenca [Hz]')
ax[0].set_ylabel('amplituda [N, m/s^2]')
ax[0].legend()
ax[1].semilogy(freq, np.abs(A/F)[segment], label='ocena FRF')
ax[1].set_xlabel('frekvenca [Hz]')
ax[1].set_ylabel('ocena FRF')
ax[1].legend();

Cenilke frekvenčne prenosne funkcije#
Pri realnih meritvah se ne moremo izogniti naključnega šuma (merilne napake) v zajetih signalih.
Ob upoštevanju odziva LTI sistemov na naključno vzbujanje z merilno napako lahko zapišemo različne cenilke frekvenčnih prenosnih funkcij:
Note
Cenilka
Note
Cenilka
Note
Kvaliteto meritve lahko ocenimo s koherenco. Ta je definirana kot:
Koherenca ima vrednosti na območju
Pozor: Izračun koherence je smiseln samo, če so spektralne gostote moči
Izračun cenilke #
Izračun spektrov: (ker bomo pri izračunu FR spektre delili med sabo, amplitudno normiranje ni nujno potrebno).
G_fa = F.conj() * A
G_ff = F.conj() * F
H1 = G_fa / G_ff
Show code cell source
plt.figure()
plt.semilogy(freq, np.abs(H1[segment]), label=f'H1, segment {segment}')
plt.semilogy(freq, np.abs(A[segment] / F[segment]), alpha=0.5, label=f'A/F, segment {segment}')
plt.xlim(0, 4000)
plt.xlabel('frekvenca [Hz]')
plt.ylabel('FRF')
plt.legend();

S povprečenjem
H1_povprečenje = np.mean(G_fa, axis=0) / np.mean(G_ff, axis=0)
Show code cell source
plt.figure()
plt.semilogy(freq, np.abs(H1_povprečenje), label=f'H1, povprečenje')
plt.semilogy(freq, np.abs(A[segment] / F[segment]), alpha=0.5, label=f'A/F, segment {segment}')
plt.xlim(0, 4000)
plt.xlabel('frekvenca [Hz]')
plt.ylabel('FRF')
plt.legend();

Izračun povprečenih spektralnih gostot moči s prekrivanjem segmentov in uporabo oken#
Welcheva metoda določi (križne in avto-) spektralne gostote moči signalov z uporabo povprečenja in oknjenja podanega signala, segmentiranega s prekrivanjem posameznih odsekov.
V Scipy je implementirana v scipy.signal.csd
(križno-spektralna gostota moči) in scipy.signal.welch
(avtospektralna gostota moči). Glavni parametri scipy.signal.csd
so:
x
: prvi signal.y
: drugi signal.fs
: frekvenca vzorčenja signalovx
iny
,window
: okno, ki naj se uporabi na posameznem segmentu. Privzeto"hann"
.nperseg
: dolžina posameznega odseka signala pri segmentiranju,noverlap
: število vzorcev, ki se med segmenti prekrivajo. Privzeto polovicanperseg
.
Cenilka z metodo Welch#
Ker imamo v vseh segmentih naključno vzbujanje, lahko zajete segmente združimo:
a = np.hstack(acc)
f = np.hstack(force)
a.shape
(256000,)
freq_welch, G_aa = signal.csd(a, a, nperseg=len(time)//4, fs=fs)
freq_welch, G_af = signal.csd(a, f, nperseg=len(time)//4, fs=fs)
(Za ustrezno normiranje križnih in avtospektrov glejte ostale argumente funkcije csd
. Pri izračunu FRF normiranje spektrov na rezultat ne vpliva.)
H2_welch = G_aa / G_af
H2_welch.shape, freq_welch.shape
((6401,), (6401,))
Show code cell source
plt.figure()
plt.semilogy(freq_welch, np.abs(H2_welch), label=f'H2, Welch')
plt.xlim(0, 4000)
plt.xlabel('frekvenca [Hz]')
plt.ylabel('FRF')
plt.legend();
