AKADEMIA
GÓRNICZO-HUTNICZA
Identyfikacja i Analiza Sygnałów
PROJEKT
Kulawik Bartosz
GR. 23, IMiR
AKADEMIA
GÓRNICZO-HUTNICZA
Identyfikacja i Analiza Sygnałów
PROJEKT
Kulawik Bartosz
GR. 23, IMiR
1. Schemat układu:
Parametry układu:
Masy: m
1
=4.5, m
2
=3, m
3
=4.5, m
4
=8, m
5
=4, m
6
=4.5
Współczynniki tłumienia: c
1
=5, c
2
=4, c
3
=1, c
4
=18, c
5
=10, c
6
=12, c
7
=4, c
8
=16,
c
9
=3
Współczynniki sprężystości: k
1
=19000, k
2
=21000, k
3
=18000, k
4
=11000,
k
5
=15000, k
6
=14000, k
7
=15000, k
8
=17000, k
9
=19000
2. Równania dynamiki w postaci równań czasowych
̈
+ ̇
+
+
(
−
̇
̇ ) +
(
−
) = 0
̈
+
(
−
̇
̇ ) +
(
−
) +
(
−
̇
̇ ) +
(
−
) = 0
̈
+ ̇
+
+
(
−
̇
̇ ) +
(
−
) =
̈
+
(
−
̇
̇ ) +
(
−
) +
(
−
̇
̇ ) +
(
−
) +
(
−
̇
̇ )
+
(
−
) +
(
−
̇
̇ ) +
(
−
) = 0
̈
+ ̇
+
+
(
−
̇
̇ ) +
(
−
) = 0
̈
+ ̇
+
+
(
−
̇
̇ ) +
(
−
) = 0
3. Równania w postaci macierzowej:
⎣
⎢
⎢
⎢
⎢
⎡
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
⎦
⎥
⎥
⎥
⎥
⎤
⎣
⎢
⎢
⎢
⎢
⎡
̈
̈
̈
̈
̈
̈ ⎦
⎥
⎥
⎥
⎥
⎤
+
⎣
⎢
⎢
⎢
⎢
⎡
+
−
0
0
0
0
−
+
0
−
0
0
0
0
+
−
0
0
0
−
−
+
+
+
−
−
0
0
0
−
+
0
0
0
0
−
0
+
⎦
⎥
⎥
⎥
⎥
⎤
⎣
⎢
⎢
⎢
⎢
⎡
̇
̇
̇
̇
̇
̇ ⎦
⎥
⎥
⎥
⎥
⎤
+
⎣
⎢
⎢
⎢
⎢
⎡
+
−
0
0
0
0
−
+
0
−
0
0
0
0
+
−
0
0
0
−
−
+
+
+
−
−
0
0
0
−
+
0
0
0
0
−
0
+
⎦
⎥
⎥
⎥
⎥
⎤
⎣
⎢
⎢
⎢
⎢
⎡
⎦
⎥
⎥
⎥
⎥
⎤
=
⎣
⎢
⎢
⎢
⎢
⎡
0
0
0
0
0 ⎦
⎥
⎥
⎥
⎥
⎤
4. Przedstawić analityczne rozwiązanie układu
Parametry układu wyznaczam rozwiązując układ równań macierzowych z
punktu 3.
Program służący do rozwiązania zagadnienia własnego w mat labie:
% Liczba stopni swobody
n=6;
% Masy w układzie
m1=4.5; m2=3; m3=4.5; m4=8;
m5=4; m6=4.5;
% Współczynniki tłumienia
c1=5; c2=4; c3=1; c4=18;
c5=10; c6=12; c7=4; c8=16; c9=3;
% Współczynniki sztywności
k1=19000; k2=21000; k3=18000; k4=11000;
k5=15000; k6=14000; k7=15000; k8=17000; k9=19000;
%macierz mas
M = [m1 0 0 0 0 0
0 m2 0 0 0 0
0 0 m3 0 0 0
0 0 0 m4 0 0
0 0 0 0 m5 0
0 0 0 0 0 m6];
%macierz współczynników tłumienia
C = [c1+c2, -c2, 0, 0, 0, 0
-c2, c2+c4, 0, -c4, 0, 0
0, 0, c3+c5, -c5, 0, 0
0, -c4, -c5, c4+c5+c6+c7, -c6, -c7
0, 0, 0, -c6, c8+c6, 0
0, 0, 0, -c7, 0, c9+c7];
%macierz współczynników sztywności
K = [k1+k2, -k2, 0, 0, 0, 0
-k2, k2+k4, 0, -k4, 0, 0
0, 0, k3+k5, -k5, 0, 0
0, -k4, -k5, k4+k5+k6+k7, -k6, -k7
0, 0, 0, -k6, k8+k6, 0
0, 0, 0, -k7, 0, k9+k7];
%budowanie macierzy
ZER = zeros(size(M));
A = [ZER,M;M,C];
B = [-M,ZER;ZER,K];
%rozwiązanie zgadnienia własnego
[PHI,LAMBDA]=eig(-B,A);
%czestotliwosci drgan wlasnych tlumionych
WD=imag(diag(LAMBDA))/2/pi;
%czestotliwosci drgan wlasnych
WW=sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2)/2/pi;
%tlumienie modalne
KSI=-real(diag(LAMBDA))./sqrt(imag(diag(LAMBDA)).^2+real(diag(LAMBDA)).^2);
disp(
'Czest. drgań Tłumienia '
)
disp(
'własnych: modalne:'
)
[WW KSI]
Rezultat działania skryptu:
Czest. drgań Tłumienia
własnych: modalne:
ans =
20.1077 0.0305
20.1077 0.0305
16.9611 0.0296
16.9611 0.0296
7.3003 0.0085
7.3003 0.0085
10.7026 0.0278
10.7026 0.0278
13.9259 0.0292
13.9259 0.0292
13.7330 0.0128
13.7330 0.0128
5. Za pomocą dowolnej metody przeprowadzić symulacje układu
Układ symuluje korzystając z równań stanu. Wymuszenie na masie 3,
rejestruje przemieszczenie masy 6. Skrypt:
%współczynniki sprężystości
k1=19000; k2=21000; k3=18000; k4=11000;
k5=15000; k6=14000; k7=15000; k8=17000; k9=19000;
%współczynniki tłumienia
c1=5; c2=4; c3=1; c4=18;
c5=10; c6=12; c7=4; c8=16; c9=3;
%masy
m1=4.5; m2=3; m3=4.5; m4=8;
m5=4; m6=4.5;
%równania stanu
A=[0 1 0 0 0 0 0 0 0 0 0 0;
-(k1+k2)/m1 -(c1+c2)/m1 k2/m1 c2/m1 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0;
k2/m2 c2/m2 -(k2+k4)/m2 -(c2+c4)/m2 0 0 k4/m2 c4/m2 0 0 0 0;
0 0 0 0 0 1 0 0 0 0 0 0;
0 0 0 0 -(k3+k5)/m3 -(c3+c5)/m3 k5/m3 c5/m3 0 0 0 0;
0 0 0 0 0 0 0 1 0 0 0 0;
0 0 k4/m4 c4/m4 k5/m4 c5/m4 -(k4+k5+k6+k7)/m4 -(c4+c5+c6+c7)/m4 k6/m4
c6/m4 k7/m4 c7/m4;
0 0 0 0 0 0 0 0 0 1 0 0;
0 0 0 0 0 0 k6/m5 c6/m5 -(k6+k8)/m5 -(c6+c8)/m5 0 0;
0 0 0 0 0 0 0 0 0 0 0 1;
0 0 0 0 0 0 k7/m6 c7/m6 0 0 -(k7+k9)/m6 -(c7+c9)/m6];
B=[0; 0; 0; 0; 0; 1/m3; 0; 0; 0; 0; 0; 0];
C=[1 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 0 0 0 0 0 1 0];
D=[0; 0; 0; 0; 0; 0];
sys=ss(A, B, C, D);
t=0:0.001:10;
u=10*rand(1, length(t));
Y_out=lsim(sys, u, t);
plot(t, Y_out(:,6),
'r'
,
'LineWidth'
,2)
title(
'Wymuszenie szumem'
);
xlabel(
'Czas [s]'
)
grid
on
figure
Y_out=impulse(sys,t);
plot(t, Y_out(:,6),
'r'
,
'LineWidth'
,2)
title(
'Wymuszenie impulsem'
);
xlabel(
'Czas [s]'
)
grid
on
Efekt działania skryptu:
Wymuszenie układu impulsem – odpowiedź masy 6
Wymuszenie układu szumem – odpowiedź masy 6
6. Identyfikacja parametryczna układu
Identyfikacje parametryczną przeprowadzam z wykorzystaniem modeli
ARMAX i ARX.
Skrypt w Matlabie:
sys=ss(A, B, C, D);
%wyznaczenie wsp. tlumienia i czest. drgan
[Wn,Z] = damp(sys);
t = 0:0.01:50;
u = 0.5*rand(1,length(t));
Y_out = lsim(sys,u,t);
%zbudowanie modelu ARMAX
model_armax = armax([Y_out(:,6) u'],[24 24 24 0]);
model_z = zpk(model_armax);
%wyznaczenie czest. drgan własnych i wsp. tlum.
%na podstawie modelu
disp(
'Model ARMAX'
)
[Wn_es, Z_es] = damp(model_z);
disp(
'Czestotliwosci drgan wlasnych'
)
disp(
'wyznaczone: wyestymowane:'
)
[Wn(:)/2/pi Wn_es(2:13)/2/pi*100]
disp(
'Wspolczynniki tłumienia'
)
disp(
'wyznaczone: wyestymowane:'
)
[Z(:) Z_es(2:13)]
%zbudowanie modelu ARX
model_arx=arx([Y_out(:,6) u'],[24 24 0]);
model_z=zpk(model_arx);
[Wn_es, Z_es] = damp(model_z);
disp(
'Model ARX'
)
[Wn_es, Z_es] = damp(model_z);
disp(
'Czestotliwosci drgan wlasnych'
)
disp(
'wyznaczone: wyestymowane:'
)
[Wn(:)/2/pi Wn_es(2:13)/2/pi*100]
disp(
'Wspolczynniki tłumienia'
)
disp(
'wyznaczone: wyestymowane:'
)
[Z(:) Z_es(2:13)]
Zestawienie wyników:
Parametry wyznaczone
Parametry estymowane
modelem ARX
Parametry estymowane
modelem ARMAX
Częstotliwość
drgań własnych
Tłumienia
modalne
Częstotliwość
drgań
własnych
Tłumienia
modalne
Częstotliwość
drgań
własnych
Tłumienia
modalne
7.3003
0.0085
7.3003
0.0085
7.3003
0.0085
7.3003
0.0085
7.3003
0.0085
7.3003
0.0085
10.7026
0.0278
10.7027
0.0278
10.7030
0.0278
10.7026
0.0278
10.7027
0.0278
10.7030
0.0278
13.7330
0.0128
13.7334
0.0128
13.7316
0.0130
13.7330
0.0128
13.7334
0.0128
13.7316
0.0130
13.9259
0.0292
13.9250
0.0294
13.9236
0.0281
13.9259
0.0292
13.9250
0.0294
13.9236
0.0281
16.9611
0.0296
16.9611
0.0296
16.9613
0.0296
16.9611
0.0296
16.9611
0.0296
16.9613
0.0296
20.1077
0.0305
20.1171
0.0304
20.1077
0.0309
20.1077
0.0305
20.1171
0.0304
20.1077
0.0309
7. Identyfikacja nieparametryczna
Skrypt:
sys=ss(A, B, C, D);
Fs=1000;
t = 0:1/Fs:10;
u = 1*rand(1,length(t));
Y_out1 = lsim(sys,u,t);
xn1=Y_out1(:,6);
Y_out=impulse(sys,t);
xn=Y_out(:,6);
figure
%periodogram
ilosc=8192;
Pxx=(abs(fft(xn, ilosc)).^2)/ilosc;
plot((0:(ilosc-1))/ilosc*Fs,10*log10(Pxx),
'black'
)
xlabel(
'Czestotliwosc [Hz]'
)
ylabel(
'Widmo mocy [dB]'
)
title(
'Periodogram'
)
xlim([0 30])
grid
on
figure
w1 = flattopwin(length(xn));
y = fft(xn.*w1);
f = (0:length(y)-1)'*Fs/length(y);
ind=find(f<=Fs/2);
f2=f(ind); am2 = 20*log10(abs(y(ind)/norm(w1)));
plot(f2,am2)
xlim([0 30])
grid
on
figure
%funkcja przejścia
nfft=ilosc;
window=hanning(nfft+1);
tfestimate(u,xn1,window,2048,nfft,Fs);
xlim([0 0.025])
figure
%funkcja koherencji
window=hanning(4096);
mscohere(xn1,u,window,3000,4096, Fs);
xlim([0 0.025])
figure
%gętość widmowa mocy
pwelch(xn, window, 2048, 4096, Fs)
xlim([0 0.05])
Wykresy:
Periodogram
Periodogram uzyskany przy Fs=100, czasie symulacji t=0:1/Fs:400 i ilości próbek
32768
Periodogram zmodyfikowany oknem flattopwin
Funkcja przejścia
Funkcja koherencji
Funkcja gęstości widmowej mocy
8. Identyfikacja układu przy pomocy transformaty falkowej
sys=ss(A, B, C, D);
Fs=200;
t = 0:1/Fs:20;
u = 10*rand(1,length(t));
%Y_out = lsim(sys,u,t);
Y_out=impulse(sys,t);
xn=Y_out(:,5)';
falka=
'cmor40-4'
;
s1=0;
F_max=Fs/2;
F_max_zadana=30;
F_min_zadana=5;
%tworzenie wektora skali
while
s1<=F_max_zadana
s1=scal2frq(F_max,falka,1/Fs);
F_max=F_max-1;
end
s1=F_min_zadana+1;
F_min=1;
while
s1>=F_min_zadana
s1=scal2frq(F_min,falka,1/Fs);
F_min=F_min+1;
end
wektor_skali=F_max:1:F_min;
figure
funkcja = cwt(xn ,wektor_skali, falka,
'plot'
) ;
figure
surf(abs(funkcja)) ;
shading
interp
view([0 90])
ylabel(
'Parametry skali'
,
'FontSize'
,12)
xlabel(
'Próbki skali'
,
'FontSize'
,12)
title(
'Skalogram zlozonego sygnalu'
)
figure
rozmiar_f=size(funkcja);
time = repmat(t,rozmiar_f(1,1),1) ;
fr = scal2frq(wektor_skali,falka,1/Fs);
freq = repmat(fr,rozmiar_f(1,2),1);
surf(time, freq', abs(funkcja)) ;
shading
interp
;
view([90 0]) ;
grid
on
title(
'Skalogram zlozonego sygnalu po przeskalowaniu'
)
ylabel(
'Czestotliwosc [Hz]'
,
'FontSize'
,12)
xlabel(
'Czas'
,
'FontSize'
,12)
figure
l=log(abs(funkcja(86,:)));
plot(t,l)
title(
'wykres logarytmu obwiedni dla 7Hz'
)
p1 = polyfit(t(1,[700:3500]),l(1,[700:3500]),1);
v=p1(1)*t+p1(2);
hold
on
plot(t, v,
'--r'
)
grid
on
p1(1)/(2*pi*7)
Na skalogramie widać częstotliwości drgań własnych układu. Najbardziej widoczne
są częstotliwość 10.7 Hz, 13.73 Hz i 13.92 Hz. Dominują również na wykresie
gęstości widmowej mocy. Grzbiety częstotliwości 10.7 Hz i 16.96 Hz są już mniej
widoczne podobnie jak na wykresie gęstości widmowej mocy. Częstotliwość 20.1
Hz na skalogramie jest niewidoczna.
Skalogram nieprzeskalowany
Wykres logarytmu obwiedni dla częstotliwości 7.3 Hz
Na podstawie współczynnika nachylenia prostej można wyznaczyć współczynnik
tłumienia.
ans =
-0.0088
Wyznaczony współczynnik nie różni się od rzeczywistego.
9. Podsumowanie
Najlepsze wyniki dała analiza parametryczna, są one najbardziej zbliżone do
wyznaczonych analitycznie. Analiza parametryczna nie daje jednak informacji na
temat mocy każdej częstotliwości. Taką informacje daje analiza nieparametryczna
– wykres gęstości widmowej mocy. Analiza nieparametryczna nie pozwala na
wyznaczenie tłumień modalnych. Najwięcej informacji na temat układu daje
transformata falkowa. Na skalogramie widać częstotliwości drgań własnych
układu, można także ocenić ich moc na podstawie wysokości grzbietu. Ponadto
można zobaczyć również przebieg częstotliwości w czasie.