test-bench-apa/matlab/motion_identification.m

101 lines
2.2 KiB
Matlab

%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
addpath('./mat/');
% Load Data
% The data from the "noise test" and the identification test are loaded.
ht = load('huddle_test.mat', 't', 'y');
load('apa95ml_5kg_Amp_E505.mat', 't', 'um', 'y');
% Any offset value is removed:
um = detrend(um, 0); % Generated DAC Voltage [V]
y = detrend(y , 0); % Mass displacement [m]
ht.y = detrend(ht.y, 0);
% Now we add a factor 10 to take into account the gain of the voltage amplifier.
Va = 10*um;
% Comparison of the PSD with Huddle Test
Ts = t(end)/(length(t)-1);
Fs = 1/Ts;
win = hanning(ceil(1*Fs));
[pxx, f] = pwelch(y, win, [], [], Fs);
[pht, ~] = pwelch(ht.y, win, [], [], Fs);
figure;
hold on;
plot(f, sqrt(pxx), 'DisplayName', '5kg');
plot(f, sqrt(pht), 'DisplayName', 'Huddle Test');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
legend('location', 'northeast');
xlim([1, Fs/2]);
% Compute TF estimate and Coherence
[tf_est, f] = tfestimate(Va, -y, win, [], [], 1/Ts);
[co_est, ~] = mscohere( Va, -y, win, [], [], 1/Ts);
figure;
hold on;
plot(f, co_est, 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([10, 5e3]);
% #+name: fig:apa95ml_5kg_PI_coh
% #+caption: Coherence
% #+RESULTS:
% [[file:figs/apa95ml_5kg_PI_coh.png]]
% Comparison with the FEM model
load('mat/fem_simscape_models.mat', 'Ghm');
freqs = logspace(0, 4, 1000);
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = nexttile([2,1]);
hold on;
plot(f, abs(tf_est), 'DisplayName', 'Identification')
plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), 'DisplayName', 'FEM')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude [m/V]'); set(gca, 'XTickLabel', []);
legend('location', 'northeast')
hold off;
ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(tf_est))
plot(freqs, 180/pi*angle(squeeze(freqresp(Ghm, freqs, 'Hz'))))
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylim([-180, 180]);
yticks(-180:90:180);
linkaxes([ax1,ax2], 'x');
xlim([10, 5e3]);