identification tests

This commit is contained in:
Operator Cad 2020-10-23 16:18:20 +02:00
parent d6e261a6bc
commit f3c087279c
5 changed files with 17 additions and 66 deletions

BIN
mat/int_enc_id_noise.mat Normal file

Binary file not shown.

Binary file not shown.

View File

@ -13,7 +13,7 @@ encoder = data(:, 2);
u = data(:, 3);
t = data(:, 4);
save('./mat/int_enc_comp.mat', 'interferometer', 'encoder', 'u' , 't');
save('./mat/int_enc_id_noise_bis.mat', 'interferometer', 'encoder', 'u' , 't');
%%
interferometer = detrend(interferometer, 0);
@ -43,26 +43,18 @@ run setup;
win = hann(ceil(10/Ts));
[tf_geo1_est, f] = tfestimate(d, geo_1, win, [], [], 1/Ts);
[co_geo1_est, ~] = mscohere(d, geo_1, win, [], [], 1/Ts);
[tf_i_est, f] = tfestimate(u, interferometer, win, [], [], 1/Ts);
[co_i_est, ~] = mscohere(u, interferometer, win, [], [], 1/Ts);
[tf_geo2_est, ~] = tfestimate(d, geo_2, win, [], [], 1/Ts);
[co_geo2_est, ~] = mscohere(d, geo_2, win, [], [], 1/Ts);
[tf_acc1_est, ~] = tfestimate(d, encoder, win, [], [], 1/Ts);
[co_acc1_est, ~] = mscohere(d, encoder, win, [], [], 1/Ts);
[tf_acc2_est, ~] = tfestimate(d, acc_2, win, [], [], 1/Ts);
[co_acc2_est, ~] = mscohere(d, acc_2, win, [], [], 1/Ts);
[tf_e_est, ~] = tfestimate(u, encoder, win, [], [], 1/Ts);
[co_e_est, ~] = mscohere(u, encoder, win, [], [], 1/Ts);
%%
figure;
hold on;
plot(f, co_geo1_est, '-')
plot(f, co_geo2_est, '-')
plot(f, co_acc1_est, '-')
plot(f, co_acc2_est, '-')
plot(f, co_i_est, '-')
plot(f, co_e_est, '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
@ -71,58 +63,18 @@ hold off;
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_geo1_est), '-', 'DisplayName', 'Geo1')
plot(f, abs(tf_i_est), '-', 'DisplayName', 'Int')
plot(f, abs(tf_e_est), '-', 'DisplayName', 'Enc')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*unwrap(angle(tf_geo1_est)), '-')
plot(f, 180/pi*unwrap(angle(tf_i_est)), '-')
plot(f, 180/pi*unwrap(angle(tf_e_est)), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
%%
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_acc1_est), '-', 'DisplayName', 'Acc1')
plot(f, abs(tf_acc2_est), '-', 'DisplayName', 'Acc2')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*unwrap(angle(tf_acc1_est)), '-')
plot(f, 180/pi*unwrap(angle(tf_acc2_est)), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
%%
win = hann(ceil(10/Ts));
[p_acc_1, f] = pwelch(encoder, win, [], [], 1/Ts);
[co_acc12, ~] = mscohere(encoder, acc_2, win, [], [], 1/Ts);
[p_geo_1, ~] = pwelch(geo_1, win, [], [], 1/Ts);
[co_geo12, ~] = mscohere(geo_1, geo_2, win, [], [], 1/Ts);
p_acc_N = p_acc_1.*(1 - co_acc12);
p_geo_N = p_geo_1.*(1 - co_geo12);
figure;
hold on;
plot(f, sqrt(p_acc_N)./abs(tf_acc1_est));
plot(f, sqrt(p_geo_N)./abs(tf_geo1_est));
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD'); xlabel('Frequency [Hz]');

13
setup.m
View File

@ -3,7 +3,7 @@ s = tf('s');
Ts = 1e-4; % [s]
%% Pre-Filter
G_pf = 1/(1 + s/2/pi/10);
G_pf = 1/(1 + s/2/pi/20);
G_pf = c2d(G_pf, Ts, 'tustin');
% %% Force Sensor Filter (HPF)
@ -16,9 +16,8 @@ G_pf = c2d(G_pf, Ts, 'tustin');
% Kiff = c2d(Kiff, Ts, 'tustin');
%
% %% Excitation Signal
% Tsim = 180; % Excitation time + Measurement time [s]
%
% t = 0:Ts:Tsim;
% % u_exc = timeseries(chirp(t, 0.1, Tsim, 1e3, 'logarithmic'), t);
% u_exc = timeseries(chirp(t, 40, Tsim, 400, 'logarithmic'), t);
% % u_exc = timeseries(y_v, t);
Tsim = 100; % Excitation time + Measurement time [s]
t = 0:Ts:Tsim;
u_exc = timeseries(chirp(t, 10, Tsim, 40, 'logarithmic'), t);
% u_exc = timeseries(y_v, t);

Binary file not shown.