phd-micro-station-modal-analysis/matlab/modal_1_meas_setup.m

%% Clear Workspace and Close figures
clear; close all; clc;

%% Intialize Laplace variable
s = zpk('s');

%% Path for functions, data and scripts
addpath('./mat/'); % Path for data

%% Colors for the figures
colors = colororder;

% Location of the Accelerometers
% <<ssec:modal_accelerometers>>

% The location of the accelerometers fixed to the micro-station is essential as it defines where the dynamics is measured.
% A total of 23 accelerometers are fixed to the six key stages of the micro station: the lower and upper granites, the translation stage, the tilt stage, the spindle and the micro hexapod.
% The position of the accelerometers are visually shown on a CAD model in Figure ref:fig:modal_location_accelerometers and their precise locations with respect to a frame located at the point of interest are summarized in Table ref:tab:modal_position_accelerometers.
% Pictures of the accelerometers fixed to the translation stage and to the micro-hexapod are shown in Figure ref:fig:modal_accelerometer_pictures.

% As all key stages of the micro-station are foreseen to behave as solid bodies, only 6 acrshort:dof can be considered per solid body.
% However, it was chosen to use four 3-axis accelerometers (i.e. 12 measured acrshort:dof) for each considered solid body to have some redundancy and to be able to verify the solid body assumption (see Section ref:ssec:modal_solid_body_assumption).

% #+attr_latex: :options [t]{0.60\linewidth}
% #+begin_minipage
% #+name: fig:modal_location_accelerometers
% #+caption: Position of the accelerometers
% #+attr_latex: :width 0.99\linewidth :float nil
% [[file:figs/modal_location_accelerometers.png]]
% #+end_minipage
% \hfill
% #+attr_latex: :options [b]{0.38\linewidth}
% #+begin_minipage
% #+begin_scriptsize
% #+name: tab:modal_position_accelerometers
% #+caption: Positions in mm
% #+attr_latex: :environment tabularx :width \linewidth :placement [b] :align Xccc
% #+attr_latex: :booktabs t :float nil :center nil
% #+RESULTS:
% |              |  $x$ |  $y$ |  $z$ |
% |--------------+------+------+------|
% | Low. Granite | -730 | -526 | -951 |
% | Low. Granite | -735 |  814 | -951 |
% | Low. Granite |  875 |  799 | -951 |
% | Low. Granite |  865 | -506 | -951 |
% | Up. Granite  | -320 | -446 | -786 |
% | Up. Granite  | -480 |  534 | -786 |
% | Up. Granite  |  450 |  534 | -786 |
% | Up. Granite  |  295 | -481 | -786 |
% | Translation  | -475 | -414 | -427 |
% | Translation  | -465 |  407 | -427 |
% | Translation  |  475 |  424 | -427 |
% | Translation  |  475 | -419 | -427 |
% | Tilt         | -385 | -300 | -417 |
% | Tilt         | -420 |  280 | -417 |
% | Tilt         |  420 |  280 | -417 |
% | Tilt         |  380 | -300 | -417 |
% | Spindle      | -155 |  -90 | -594 |
% | Spindle      |    0 |  180 | -594 |
% | Spindle      |  155 |  -90 | -594 |
% | Hexapod      |  -64 |  -64 | -270 |
% | Hexapod      |  -64 |   64 | -270 |
% | Hexapod      |   64 |   64 | -270 |
% | Hexapod      |   64 |  -64 | -270 |
% #+end_scriptsize
% #+end_minipage

% #+name: fig:modal_accelerometer_pictures
% #+caption: Accelerometers fixed on the micro-station stages
% #+attr_latex: :options [htbp]
% #+begin_figure
% #+attr_latex: :caption \subcaption{\label{fig:modal_accelerometers_ty} $T_y$ stage}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :height 6cm
% [[file:figs/modal_accelerometers_ty.jpg]]
% #+end_subfigure
% #+attr_latex: :caption \subcaption{\label{fig:modal_accelerometers_hexapod} Micro-Hexapod}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :height 6cm
% [[file:figs/modal_accelerometers_hexapod.jpg]]
% #+end_subfigure
% #+end_figure


%% Load Accelerometer positions
acc_pos = readtable('mat/acc_pos.txt', 'ReadVariableNames', false);
acc_pos = table2array(acc_pos(:, 1:4));
[~, i] = sort(acc_pos(:, 1));
acc_pos = acc_pos(i, 2:4);

% Force and Response signals
% <<ssec:modal_measured_signals>>

% The force sensor of the instrumented hammer and the accelerometers signals are shown in the time domain in Figure ref:fig:modal_raw_meas.
% Sharp "impacts" can be seen for the force sensor, indicating wide frequency band excitation.
% For the accelerometer, a much more complex signal can be observed, indicating complex dynamics.

% The "normalized" acrfull:asd of the two signals are computed and shown in Figure ref:fig:modal_asd_acc_force.
% Conclusions based on the time domain signals can be clearly seen in the frequency domain (wide frequency content for the force signal and complex dynamics for the accelerometer).
% Similar results are obtained for all the measured frequency response functions.


%% Load raw data
meas1_raw = load('mat/meas_raw_1.mat');

% Sampling Frequency [Hz]
Fs = 1/meas1_raw.Track1_X_Resolution;

% Time just before the impact occurs [s]
impacts = [5.937, 11.228, 16.681, 22.205, 27.350, 32.714, 38.115, 43.888, 50.407]-0.01;

% Time vector [s]
time = linspace(0, meas1_raw.Track1_X_Resolution*length(meas1_raw.Track1), length(meas1_raw.Track1));

%% Raw measurement of the Accelerometer
figure;
hold on;
plot(time-22.2, meas1_raw.Track2, 'DisplayName', '$X_{j}$ [$m/s^2$]');
plot(time-22.2, 1e-3*meas1_raw.Track1, 'DisplayName', '$F_{k}$ [kN]');
hold off;
xlabel('Time [s]');
ylabel('Amplitude');
xlim([0, 0.2]);
ylim([-2, 2]);
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);

%% Frequency Analysis
Nfft = floor(5.0*Fs); % Number of frequency points
win = hanning(Nfft); % Windowing
Noverlap = floor(Nfft/2); % Overlap for frequency analysis

%% Comnpute the power spectral density of the force and acceleration
[pxx_force, f] = pwelch(meas1_raw.Track1, win, Noverlap, Nfft, Fs);
[pxx_acc,   ~] = pwelch(meas1_raw.Track2, win, Noverlap, Nfft, Fs);

%% Normalized Amplitude Spectral Density of the measured force and acceleration
figure;
hold on;
plot(f, sqrt(pxx_acc./max(pxx_acc(f<200))), 'DisplayName', '$X_{j}$');
plot(f, sqrt(pxx_force./max(pxx_force(f<200))), 'DisplayName', '$F_{k}$');
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Normalized Spectral Density');
xlim([0, 200]);
xticks([0:20:200]);
ylim([0, 1])
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);



% #+name: fig:modal_raw_meas_asd
% #+caption: Raw measurement of the acceleromter (blue) and of the force sensor at the Hammer tip (red) (\subref{fig:modal_raw_meas}). Computed Amplitude Spectral Density of the two signals (normalized) (\subref{fig:modal_asd_acc_force})
% #+attr_latex: :options [htbp]
% #+begin_figure
% #+attr_latex: :caption \subcaption{\label{fig:modal_raw_meas}Time domain signals}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :width 0.95\linewidth
% [[file:figs/modal_raw_meas.png]]
% #+end_subfigure
% #+attr_latex: :caption \subcaption{\label{fig:modal_asd_acc_force}Amplitude Spectral Density (normalized)}
% #+attr_latex: :options {0.49\textwidth}
% #+begin_subfigure
% #+attr_latex: :width 0.95\linewidth
% [[file:figs/modal_asd_acc_force.png]]
% #+end_subfigure
% #+end_figure

% The frequency response function $H_{jk}$ from the applied force $F_{k}$ to the measured acceleration $X_j$ is then computed and shown Figure ref:fig:modal_frf_acc_force.
% The quality of the obtained data can be estimated using the /coherence/ function, which is shown in Figure ref:fig:modal_coh_acc_force.
% Good coherence is obtained from $20\,\text{Hz}$ to $200\,\text{Hz}$ which corresponds to the frequency range of interest.


%% Compute the transfer function and Coherence
[G1,   f] = tfestimate(meas1_raw.Track1, meas1_raw.Track2, win, Noverlap, Nfft, Fs);
[coh1, ~] = mscohere(  meas1_raw.Track1, meas1_raw.Track2, win, Noverlap, Nfft, Fs);

%% Frequency Response Function between the force and the acceleration
figure;
plot(f, abs(G1));
xlabel('Frequency [Hz]'); ylabel('FRF [$m/s^2/N$]')
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'log');
xlim([0, 200]);
xticks([0:20:200]);

%% Frequency Response Function between the force and the acceleration
figure;
plot(f, coh1);
xlabel('Frequency [Hz]'); ylabel('Coherence [-]')
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlim([0, 200]); ylim([0,1]);
xticks([0:20:200]);