Update analysis => FRF in the global frame

modal-analysis/frf_processing.org

@@ -0,0 +1,680 @@
+#+TITLE: Modal Analysis - Processing of FRF
+:DRAWER:
+#+STARTUP: overview
+
+#+LANGUAGE: en
+#+EMAIL: dehaeze.thomas@gmail.com
+#+AUTHOR: Dehaeze Thomas
+
+#+HTML_LINK_HOME: ../index.html
+#+HTML_LINK_UP: ./index.html
+
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
+#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
+#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
+#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
+#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
+#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>
+
+#+HTML_MATHJAX: align: center tagside: right font: TeX
+
+#+PROPERTY: header-args:matlab  :session *MATLAB*
+#+PROPERTY: header-args:matlab+ :comments org
+#+PROPERTY: header-args:matlab+ :results none
+#+PROPERTY: header-args:matlab+ :exports both
+#+PROPERTY: header-args:matlab+ :eval no-export
+#+PROPERTY: header-args:matlab+ :output-dir figs
+
+#+PROPERTY: header-args:shell  :eval no-export
+
+#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
+#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
+#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
+#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
+#+PROPERTY: header-args:latex+ :results raw replace :buffer no
+#+PROPERTY: header-args:latex+ :eval no-export
+#+PROPERTY: header-args:latex+ :exports both
+#+PROPERTY: header-args:latex+ :mkdirp yes
+#+PROPERTY: header-args:latex+ :output-dir figs
+:END:
+
+The measurements have been conducted and we have computed the $n \times p \times q$ Frequency Response Functions Matrix with:
+- $n$: the number of measurements: $23 \times 3 = 69$ (23 accelerometers measuring 3 directions each)
+- $p$: the number of excitation inputs: $3$
+- $q$: the number of frequency points $\omega_i$
+
+However, in our model, we only consider 6 solid bodies, namely:
+- Bottom Granite
+- Top Granite
+- Translation Stage
+- Tilt Stage
+- Spindle
+- Hexapod
+
+Thus, we are only interested in $6 \times 6 = 36$ degrees of freedom.
+
+We here process the FRF matrix to go from the 69 measured DOFs to the wanted 36 DOFs.
+
+* Matlab Init                                               :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+* Importation of measured FRF curves
+We load the measured FRF and Coherence matrices.
+We also load the geometric parameters of the station: solid bodies considered and the position of the accelerometers.
+
+#+begin_src matlab
+  load('./mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
+  load('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
+#+end_src
+
+* From accelerometer DOFs to solid body DOFs - Mathematics
+Let's consider the schematic shown on figure [[fig:local_to_global_coordinates]] where we are measuring the motion of a (supposed) solid body at 4 distinct points in x-y-z.
+
+The goal here is to link these $4 \times 3 = 12$ measurements to the 6 DOFs of the solid body expressed in the frame $\{O\}$.
+
+#+begin_src latex :file local_to_global_coordinates.pdf :post pdf2svg(file=*this*, ext="png") :exports results
+  \newcommand\irregularcircle[2]{% radius, irregularity
+    \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
+    +(0:\len pt)
+    \foreach \a in {10,20,...,350}{
+      \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
+      -- +(\a:\len pt)
+    } -- cycle
+  }
+  \begin{tikzpicture}
+    \draw[rounded corners=1mm, fill=blue!30!white] (0, 0) \irregularcircle{3cm}{1mm};
+
+    \node[] (origin) at (4, -1) {$\bullet$};
+    \begin{scope}[shift={(origin)}]
+      \def\axissize{0.8cm}
+      \draw[->] (0, 0) -- ++(\axissize, 0) node[above left]{$x$};
+      \draw[->] (0, 0) -- ++(0, \axissize) node[below right]{$y$};
+      \draw[fill, color=black] (0, 0) circle (0.05*\axissize);
+      \node[draw, circle, inner sep=0pt, minimum size=0.4*\axissize, label=left:$z$] (yaxis) at (0, 0){};
+      \node[below right] at (0, 0){$\{O\}$};
+    \end{scope}
+
+    \coordinate[] (p1) at (-1.5, -1.5);
+    \coordinate[] (p2) at (-1.5,  1.5);
+    \coordinate[] (p3) at ( 1.5,  1.5);
+    \coordinate[] (p4) at ( 1.5, -1.5);
+
+    \draw[->] (p1)node[]{$\bullet$}node[above]{$p_1$} -- ++(1, 0.5)node[right]{$v_1$};
+    \draw[->] (p2)node[]{$\bullet$}node[above]{$p_2$} -- ++(-0.5, 1)node[right]{$v_2$};
+    \draw[->] (p3)node[]{$\bullet$}node[above]{$p_3$} -- ++(1, 0.5)node[right]{$v_3$};
+    \draw[->] (p4)node[]{$\bullet$}node[above]{$p_4$} -- ++(0.5, 1)node[right]{$v_4$};
+  \end{tikzpicture}
+#+end_src
+
+#+name: fig:local_to_global_coordinates
+#+caption: Schematic of the measured motions of a solid body
+#+RESULTS:
+[[file:figs/local_to_global_coordinates.png]]
+
+From the figure [[fig:local_to_global_coordinates]], we can write:
+\begin{align*}
+  \vec{v}_1 &= \vec{v} + \Omega \vec{p}_1\\
+  \vec{v}_2 &= \vec{v} + \Omega \vec{p}_2\\
+  \vec{v}_3 &= \vec{v} + \Omega \vec{p}_3\\
+  \vec{v}_4 &= \vec{v} + \Omega \vec{p}_4
+\end{align*}
+
+With
+\begin{equation}
+  \Omega = \begin{bmatrix}
+  0 & -\Omega_z & \Omega_y \\
+  \Omega_z & 0 & -\Omega_x \\
+  -\Omega_y & \Omega_x & 0
+\end{bmatrix}
+\end{equation}
+$\vec{v}$ and $\Omega$ represent to velocity and rotation of the solid expressed in the frame $\{O\}$.
+
+We can rearrange the equations in a matrix form:
+
+\begin{equation}
+\left[\begin{array}{ccc|ccc}
+  1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\
+  0 & 1 & 0 & -p_{1z} & 0 & p_{1x} \\
+  0 & 0 & 1 & p_{1y} & -p_{1x} & 0 \\ \hline
+    & \vdots &   &    & \vdots &   \\ \hline
+  1 & 0 & 0 & 0 & p_{4z} & -p_{4y} \\
+  0 & 1 & 0 & -p_{4z} & 0 & p_{4x} \\
+  0 & 0 & 1 & p_{4y} & -p_{4x} & 0
+\end{array}\right] \begin{bmatrix}
+  v_x \\ v_y \\ v_z \\ \hline \Omega_x \\ \Omega_y \\ \Omega_z
+\end{bmatrix} = \begin{bmatrix}
+  v_{1x} \\ v_{1y} \\ v_{1z} \\\hline \vdots \\\hline v_{4x} \\ v_{4y} \\ v_{4z}
+\end{bmatrix}
+\end{equation}
+
+and then we obtain the velocity and rotation of the solid in the wanted frame $\{O\}$:
+\begin{equation}
+\begin{bmatrix}
+  v_x \\ v_y \\ v_z \\ \hline \Omega_x \\ \Omega_y \\ \Omega_z
+\end{bmatrix} =
+\left[\begin{array}{ccc|ccc}
+  1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\
+  0 & 1 & 0 & -p_{1z} & 0 & p_{1x} \\
+  0 & 0 & 1 & p_{1y} & -p_{1x} & 0 \\ \hline
+    & \vdots &   &    & \vdots &   \\ \hline
+  1 & 0 & 0 & 0 & p_{4z} & -p_{4y} \\
+  0 & 1 & 0 & -p_{4z} & 0 & p_{4x} \\
+  0 & 0 & 1 & p_{4y} & -p_{4x} & 0
+\end{array}\right]^{-1} \begin{bmatrix}
+  v_{1x} \\ v_{1y} \\ v_{1z} \\\hline \vdots \\\hline v_{4x} \\ v_{4y} \\ v_{4z}
+\end{bmatrix}
+\end{equation}
+
+This inversion is equivalent to resolving a mean square problem.
+
+* What reference frame to choose?
+The question we wish here to answer is how to choose the reference frame $\{O\}$ in which the DOFs of the solid bodies are defined.
+
+The goal is to compare the motion of each solid body to see which relative DOFs between solid bodies can be neglected, that is to say, which joint between solid bodies can be regarded as perfect (and this in all the frequency range of interest).
+Ideally, we would like to have the same number of degrees of freedom than the number of identified modes.
+
+The possibles choices are:
+- *One frame for each solid body* which is located at its center of mass
+- *One common frame*, for instance located at the point of interest ($270mm$ above the Hexapod)
+- *Base located at the joint position*: this is where we want to see the motion and estimate stiffness
+
+#+name: tab:frame_comparison
+#+caption: Advantages and disadvantages for the choice of reference frame
+| Chosen Frame             | Advantages                                          | Disadvantages                                        |
+|--------------------------+-----------------------------------------------------+------------------------------------------------------|
+| Frames at CoM            | Physically, it makes more sense                     | How to compare the motion of the solid bodies?       |
+| Common Frame             | We can compare the motion of each solid body        | Small $\theta_{x, y}$ may result in large $T_{x, y}$ |
+| Frames at joint position | Directly gives which joint direction can be blocked | How to choose the joint position?                    |
+
+As the easiest choice is to choose a common frame, we start with that solution.
+
+* From accelerometer DOFs to solid body DOFs - Matlab Implementation
+First, we initialize a new FRF matrix =FRFs_O= which is an $n \times p \times q$ with:
+- $n$ is the number of DOFs of the considered 6 solid-bodies: $6 \times 6 = 36$
+- $p$ is the number of excitation inputs: $3$
+- $q$ is the number of frequency points $\omega_i$
+
+#+begin_src matlab
+  FRFs_O = zeros(length(solid_names)*6, 3, 801);
+#+end_src
+
+Then, as we know the positions of the accelerometers on each solid body, and we have the response of those accelerometers, we can use the equations derived in the previous section to determine the response of each solid body expressed in the frame $\{O\}$.
+#+begin_src matlab
+  for solid_i = 1:length(solid_names)
+    solids_i = solids.(solid_names{solid_i});
+
+    A = zeros(3*length(solids_i), 6);
+    for i = 1:length(solids_i)
+      acc_i = solids_i(i);
+
+      A(3*(i-1)+1:3*i, 1:3) = eye(3);
+      A(3*(i-1)+1:3*i, 4:6) = [ 0                  acc_pos(acc_i, 3) -acc_pos(acc_i, 2) ;
+                               -acc_pos(acc_i, 3)  0                  acc_pos(acc_i, 1) ;
+                                acc_pos(acc_i, 2) -acc_pos(acc_i, 1)  0];
+    end
+
+    for exc_dir = 1:3
+      FRFs_O((solid_i-1)*6+1:solid_i*6, exc_dir, :) = A\squeeze(FRFs((solids_i(1)-1)*3+1:solids_i(end)*3, exc_dir, :));
+    end
+  end
+#+end_src
+
+* Analysis of some FRF in the global coordinates
+First, we can compare the motions of the 6 solid bodies in one direction (figure [[fig:frf_all_bodies_one_direction]])
+
+We can also compare all the DOFs of one solid body (figure [[fig:frf_one_body_all_directions]]).
+
+#+begin_src matlab :exports none
+  exc_names = {'$F_x$', '$F_y$', '$F_z$'};
+  DOFs = {'$T_x$', '$T_y$', '$T_z$', '$\theta_x$', '$\theta_y$', '$\theta_z$'}
+  solids_i = 1:6;
+  dir_i = 1;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subaxis(2, 1, 1);
+  hold on;
+  for solid_i = solids_i
+    plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 'DisplayName', solid_names{solid_i});
+  end
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+  legend('Location', 'northwest');
+  title(sprintf('FRF between %s and %s', exc_names{exc_dir}, DOFs{dir_i}));
+
+  ax2 = subaxis(2, 1, 2);
+  hold on;
+  for solid_i = solids_i
+    plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 360)-180);
+  end
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/frf_all_bodies_one_direction.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:frf_all_bodies_one_direction
+#+CAPTION: FRFs of all the 6 solid bodies in one direction
+[[file:figs/frf_all_bodies_one_direction.png]]
+
+#+begin_src matlab :exports none
+  DOFs = {'$T_x$', '$T_y$', '$T_z$', '$\theta_x$', '$\theta_y$', '$\theta_z$'}
+  solid_i = 3;
+  dirs_i = 1:6;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for dir_i = dirs_i
+    plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 'DisplayName', DOFs{dir_i});
+  end
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+  legend('Location', 'northwest');
+  title(sprintf('Motion of %s due to %s', solid_names{solid_i}, exc_names{exc_dir}));
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for dir_i = dirs_i
+    plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 360)-180);
+  end
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/frf_one_body_all_directions.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:frf_one_body_all_directions
+#+CAPTION: FRFs of one solid body in all its DOFs
+[[file:figs/frf_one_body_all_directions.png]]
+
+* TODO How to compare the relative motion of solid bodies
+We have some of elements of the full FRF matrix:
+\[ \frac{D_{1x}}{F_x},\ \frac{D_{1y}}{F_x},\ \frac{D_{1z}}{F_x},\ \frac{D_{2x}}{F_x},\ \dots \]
+
+\[ \frac{D_{1x}}{D_{2x}} = \frac{\frac{D_{1x}}{F_x}}{\frac{D_{2x}}{F_x}} \]
+Then, if $\left| \frac{D_{1x}}{D_{2x}} \right| \approx 1$ in all the frequency band of interest, we can block the $x$ motion between the solids 1 and 2.
+
+\[ \frac{D_{2x} - D_{1x}}{D_{1x} + D_{2x}} = \frac{\frac{D_{2x}}{F_x} - \frac{D_{1x}}{F_x}}{\frac{D_{1x}}{F_x} + \frac{D_{2x}}{F_x}} \]
+
+Then if $\left| \frac{D_{2x} - D_{1x}}{D_{1x} + D_{2x}} \right| \ll 1$ in all the frequency band of interest, we can block the $x$ motion between the solids 1 and 2.
+
+* Relative Motion in the global coordinates
+Below we plot the normalized relative motion between each stage:
+\[ 0 < \frac{\left| D_{ix} - D_{jx} \right|}{|D_{ix}| + |D_{jx}|} < 1 \]
+
+#+begin_src matlab
+  DOFs = {'$T_x$', '$T_y$', '$T_z$', '$\theta_x$', '$\theta_y$', '$\theta_z$'}
+
+  dirs_i = 1:6;
+  exc_dir = 1;
+
+  figure;
+
+  for i = 2:6
+    subaxis(3, 2, i);
+    hold on;
+    for dir_i = dirs_i
+      H = (squeeze(FRFs_O((i-1)*6+dir_i, exc_dir, :))-squeeze(FRFs_O((i-2)*6+dir_i, exc_dir, :)))./(abs(squeeze(FRFs_O((i-1)*6+dir_i, exc_dir, :)))+abs(squeeze(FRFs_O((i-2)*6+dir_i, exc_dir, :))));
+      plot(freqs, abs(H));
+    end
+    hold off;
+    set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+    xlim([1, 200]); ylim([0, 1]);
+    % xlabel('Frequency [Hz]'); ylabel('Relative Motion');
+    title(sprintf('Normalized motion %s - %s', solid_names{i-1}, solid_names{i}));
+    if i > 4
+      xlabel('Frequency [Hz]');
+    else
+      set(gca, 'XTickLabel',[]);
+    end
+  end
+
+  for i = 1:length(dirs_i)
+    legend_names{i} = DOFs{dirs_i(i)};
+  end
+  lgd = legend(legend_names);
+
+  hL = subplot(3, 2, 1);
+  poshL = get(hL,'position');
+
+  set(lgd,'position', poshL);
+  axis(hL, 'off');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/relative_motion_comparison.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:relative_motion_comparison
+#+CAPTION: Relative motion between each stage
+[[file:figs/relative_motion_comparison.png]]
+
+* TODO Compare original FRF measurements to transformed FRF in the global frame
+We wish here to compare the FRF in order to verify if there is any mistake.
+
+#+begin_src matlab
+  dir_names = {'X', 'Y', 'Z', '$\theta_X$', '$\theta_Y$', '$\theta_Z$'};
+
+  solid_i = 6;
+  acc_dir_O = 1;
+  acc_dir = 1;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for i = solids.(solid_names{solid_i})
+    plot(freqs, abs(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))));
+  end
+  plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), '-k');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+  title(sprintf('%s motion measured by the Acc. vs %s motion computed in the common frame - %s', dir_names{acc_dir}, dir_names{acc_dir_O}, solid_names{solid_i}));
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for i = solids.(solid_names{solid_i})
+    plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))), 360)-180);
+  end
+  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+acc_dir_O, exc_dir, :))), 360)-180, '-k');
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+* Verify that we find the original FRF from the FRF in the global coordinates
+We have computed the Frequency Response Functions Matrix =FRFs_O= representing the response of the 6 solid bodies in their 6 DOFs.
+
+From the response of one body in its 6 DOFs, we should be able to compute the FRF of each of its accelerometer fixed to it during the measurement.
+
+We can then compare the result with the original measurements.
+This will help us to determine if:
+- the previous inversion used is correct
+- the solid body assumption is correct in the frequency band of interest
+
+#+begin_src matlab
+  FRF_recovered = zeros(size(FRFs));
+
+  % For each excitation direction
+  for exc_dir = 1:3
+    % For each solid
+    for solid_i = 1:length(solid_names)
+      v0 = squeeze(FRFs_O((solid_i-1)*6+1:(solid_i-1)*6+3, exc_dir, :));
+      W0 = squeeze(FRFs_O((solid_i-1)*6+4:(solid_i-1)*6+6, exc_dir, :));
+
+      % For each accelerometer attached to the current solid
+      for acc_i = solids.(solid_names{solid_i})
+        % We get the position of the accelerometer expressed in frame O
+        pos = acc_pos(acc_i, :)';
+        posX = [0 pos(3) -pos(2); -pos(3) 0 pos(1) ; pos(2) -pos(1) 0];
+               [0 acc_pos(i, 3) -acc_pos(i, 2) ; -acc_pos(i, 3) 0 acc_pos(i, 1) ; acc_pos(i, 2) -acc_pos(i, 1) 0]
+
+        FRF_recovered(3*(acc_i-1)+1:3*(acc_i-1)+3, exc_dir, :) = v0 + posX*W0;
+      end
+    end
+  end
+#+end_src
+
+We then compare the original FRF measured for each accelerometer with the recovered FRF from the global FRF matrix in the common frame.
+
+The FRF for the 4 accelerometers on the Hexapod are compared on figure [[fig:recovered_frf_comparison_hexa]].
+All the FRF are matching very well in all the frequency range displayed.
+
+The FRF for accelerometers located on the translation stage are compared on figure [[fig:recovered_frf_comparison_ty]].
+The FRF are matching well until 100Hz.
+
+#+begin_src matlab :exports none
+  exc_names = {'$F_x$', '$F_y$', '$F_z$'};
+  DOFs = {'$T_x$', '$T_y$', '$T_z$', '$\theta_x$', '$\theta_y$', '$\theta_z$'}
+
+  solid_i = 6;
+  exc_dir = 1;
+
+  accs_i = solids.(solid_names{solid_i});
+
+  figure;
+
+  for i = 1:length(accs_i)
+    acc_i = accs_i(i);
+
+    subaxis(2, 2, i);
+
+    hold on;
+    for dir_i = 1:3
+      plot(freqs, abs(squeeze(FRFs(3*(acc_i-1)+dir_i, exc_dir, :))), '-', 'DisplayName', DOFs{dir_i});
+    end
+    set(gca,'ColorOrderIndex',1)
+    for dir_i = 1:3
+      plot(freqs, abs(squeeze(FRF_recovered(3*(acc_i-1)+dir_i, exc_dir, :))), '--', 'HandleVisibility', 'off');
+    end
+    hold off;
+
+    set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+    if i > 2
+      xlabel('Frequency [Hz]');
+    else
+      set(gca, 'XTickLabel',[]);
+    end
+
+    if rem(i, 2) == 1
+      ylabel('Amplitude');
+    end
+
+    xlim([1, 200]);
+    title(sprintf('Accelerometer %i', accs_i(i)));
+    legend('location', 'northwest');
+  end
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/recovered_frf_comparison_hexa.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:recovered_frf_comparison_hexa
+#+CAPTION: Comparison of the original FRF with the recovered ones - Hexapod
+[[file:figs/recovered_frf_comparison_hexa.png]]
+
+#+begin_src matlab :exports none
+  exc_names = {'$F_x$', '$F_y$', '$F_z$'};
+  DOFs = {'$T_x$', '$T_y$', '$T_z$', '$\theta_x$', '$\theta_y$', '$\theta_z$'}
+
+  solid_i = 3;
+  exc_dir = 1;
+
+  accs_i = solids.(solid_names{solid_i});
+
+  figure;
+
+  for i = 1:length(accs_i)
+    acc_i = accs_i(i);
+
+    subaxis(2, 2, i);
+
+    hold on;
+    for dir_i = 1:3
+      plot(freqs, abs(squeeze(FRFs(3*(acc_i-1)+dir_i, exc_dir, :))), '-', 'DisplayName', DOFs{dir_i});
+    end
+    set(gca,'ColorOrderIndex',1)
+    for dir_i = 1:3
+      plot(freqs, abs(squeeze(FRF_recovered(3*(acc_i-1)+dir_i, exc_dir, :))), '--', 'HandleVisibility', 'off');
+    end
+    hold off;
+
+    set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+    if i > 2
+      xlabel('Frequency [Hz]');
+    else
+      set(gca, 'XTickLabel',[]);
+    end
+
+    if rem(i, 2) == 1
+      ylabel('Amplitude');
+    end
+
+    xlim([1, 200]);
+    title(sprintf('Accelerometer %i', accs_i(i)));
+    legend('location', 'northwest');
+  end
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/recovered_frf_comparison_ty.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:recovered_frf_comparison_ty
+#+CAPTION: Comparison of the original FRF with the recovered ones - Ty
+[[file:figs/recovered_frf_comparison_ty.png]]
+
+
+#+begin_important
+  The reduction of the number of degrees of freedom from 69 (23 accelerometers with each 3DOF) to 36 (6 solid bodies with 6 DOF) seems to work well.
+
+  This confirms the fact that the stages are indeed behaving as a solid body in the frequency band of interest.
+  This valid the fact that a multi-body model can be used to represent the dynamics of the micro-station.
+#+end_important
+
+* Importation of measured FRF curves                        :noexport:ignore:
+There are 24 measurements files corresponding to 24 series of impacts:
+- 3 directions, 8 sets of 3 accelerometers
+For each measurement file, the FRF and coherence between the impact and the 9 accelerations measured.
+
+In reality: 4 sets of 10 things
+
+#+begin_src matlab
+  a = load('mat/meas_frf_coh_1.mat');
+#+end_src
+
+#+begin_src matlab
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Mod)
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+  title(sprintf('From %s, to %s', FFT1_AvXSpc_2_1_RfName, FFT1_AvXSpc_2_1_RpName))
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Phas)
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+* Analysis of some FRFs                                     :noexport:ignore:
+#+begin_src matlab
+  acc_i = 3;
+  acc_dir = 1;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(freqs, abs(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))), 360)-180);
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+
+#+begin_src matlab
+  figure;
+  hold on;
+  for i = 1:3*n_acc
+    plot(freqs, squeeze(COHs(i, 1, :)), 'color', [0, 0, 0, 0.2]);
+  end
+  hold off;
+  xlabel('Frequency [Hz]');
+  ylabel('Coherence [\%]');
+#+end_src
+
+
+Composite Response Function.
+
+We here sum the norm instead of the complex numbers.
+
+#+begin_src matlab
+  HHx = squeeze(sum(abs(FRFs(:, 1, :))));
+  HHy = squeeze(sum(abs(FRFs(:, 2, :))));
+  HHz = squeeze(sum(abs(FRFs(:, 3, :))));
+  HH  = squeeze(sum([HHx, HHy, HHz], 2));
+#+end_src
+
+#+begin_src matlab
+  exc_dir = 3;
+
+  figure;
+  hold on;
+  for i = 1:3*n_acc
+    plot(freqs, abs(squeeze(FRFs(i, exc_dir, :))), 'color', [0, 0, 0, 0.2]);
+  end
+  plot(freqs, abs(HHx));
+  plot(freqs, abs(HHy));
+  plot(freqs, abs(HHz));
+  plot(freqs, abs(HH), 'k');
+  hold off;
+  set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude');
+  xlim([1, 200]);
+#+end_src

modal-analysis/mat/acc_pos.txt

@@ -17,7 +17,7 @@
      7  4.2000e-001  2.8000e-001 -4.1680e-001     
      6 -4.2000e-001  2.8000e-001 -4.1680e-001     
      5 -3.8500e-001 -3.0000e-001 -4.1680e-001     
-     4  6.4000e-002 -6.4000e-002 -2.9600e-001     
-     3  6.4000e-002  6.4000e-002 -2.9600e-001     
-     2 -6.4000e-002  6.4000e-002 -2.9600e-001     
-     1 -6.4000e-002 -6.4000e-002 -2.9600e-001     
+     4  6.4000e-002 -6.4000e-002 -2.7000e-001     
+     3  6.4000e-002  6.4000e-002 -2.7000e-001     
+     2 -6.4000e-002  6.4000e-002 -2.7000e-001     
+     1 -6.4000e-002 -6.4000e-002 -2.7000e-001     

modal-analysis/mathematical_model.org

@@ -60,3 +60,35 @@ However, each of the DOF of the system may not be relevant for the modes present
 For instance, the translation stage may not vibrate in the Z direction for all the modes identified. Then, we can block this DOF and this simplifies the model.
 
 The modal identification done here will thus permit us to determine *which DOF can be neglected*.
+
+* Some notes about constraining the number of degrees of freedom
+We want to have the two eigen matrices.
+
+They should have the same size $n \times n$ where $n$ is the number of modes as well as the number of degrees of freedom.
+Thus, if we consider 21 modes, we should restrict our system to have only 21 DOFs.
+
+Actually, we are measured 6 DOFs of 6 solids, thus we have 36 DOFs.
+
+From the mode shapes animations, it seems that in the frequency range of interest, the two marbles can be considered as one solid.
+We thus have 5 solids and 30 DOFs.
+
+
+In order to determine which DOF can be neglected, two solutions seems possible:
+- compare the mode shapes
+- compare the FRFs
+
+The question is: in which base (frame) should be express the modes shapes and FRFs?
+Is it meaningful to compare mode shapes as they give no information about the amplitudes of vibration?
+
+
+| Stage   | Motion DOFs | Parasitic DOF | Total DOF | Description of DOF |
+|---------+-------------+---------------+-----------+--------------------|
+| Granite |           0 |             3 |         3 |                    |
+| Ty      |           1 |             2 |         3 | Ty, Rz             |
+| Ry      |           1 |             2 |         3 | Ry,                |
+| Rz      |           1 |             2 |         3 | Rz, Rx, Ry         |
+| Hexapod |           6 |             0 |         6 | Txyz, Rxyz         |
+|---------+-------------+---------------+-----------+--------------------|
+|         |           9 |             9 |        18 |                    |
+#+TBLFM: $4=vsum($2..$3)
+#+TBLFM: @>$2..$>=vsum(@I..@II)

modal-analysis/matlab/modal_frf_coh.m

@@ -4,6 +4,21 @@ clear; close all; clc;
 %% Intialize Laplace variable
 s = zpk('s');
 
+
+
+% We then import that on =matlab=, and sort them.
+
+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);
+
+
+
+% The positions of the sensors relative to the point of interest are shown below.
+
+data2orgtable([[1:23]', 1000*acc_pos], {}, {'ID', 'x [mm]', 'y [mm]', 'z [mm]'}, ' %.0f ');
+
 % Windowing
 % Windowing is used on the force and response signals.
 
@@ -301,9 +316,6 @@ for i = 1:n_meas
       meas_factor = meas_factor*(-1);
     end
 
-    % FRFs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
-    % COHs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
-
     FRFs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
     COHs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
   end
@@ -315,3 +327,30 @@ freqs = meas.FFT1_Coh_10_1_RMS_X_Val;
 % And we save the obtained FRF matrix and Coherence matrix in a =.mat= file.
 
 save('./mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
+
+% Solid Bodies considered for further analysis
+% We consider the following solid bodies for further analysis:
+% - Bottom Granite
+% - Top Granite
+% - Translation Stage
+% - Tilt Stage
+% - Spindle
+% - Hexapod
+
+% We create a =matlab= structure =solids= that contains the accelerometers ID connected to each solid bodies (as shown on figure [[fig:nass-modal-test]]).
+
+solids = {};
+solids.granite_bot = [17, 18, 19, 20];
+solids.granite_top = [13, 14, 15, 16];
+solids.ty          = [9, 10, 11, 12];
+solids.ry          = [5, 6, 7, 8];
+solids.rz          = [21, 22, 23];
+solids.hexa        = [1, 2, 3, 4];
+
+solid_names = fields(solids);
+
+
+
+% Finally, we save that into a =.mat= file.
+
+save('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');

modal-analysis/measurement.org

@@ -29,6 +29,16 @@
 #+PROPERTY: header-args:matlab+ :mkdirp yes
 
 #+PROPERTY: header-args:shell  :eval no-export
+
+#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
+#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
+#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
+#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
+#+PROPERTY: header-args:latex+ :results raw replace :buffer no
+#+PROPERTY: header-args:latex+ :eval no-export
+#+PROPERTY: header-args:latex+ :exports both
+#+PROPERTY: header-args:latex+ :mkdirp yes
+#+PROPERTY: header-args:latex+ :output-dir figs
 :END:
 
 * ZIP file containing the data and matlab files                      :ignore:
@@ -57,6 +67,9 @@
   <<matlab-init>>
 #+end_src
 
+* Goal
+The goal is to measure the dynamic of the Micro-Station and to extract Frequency Response Functions.
+
 * Instrumentation Used
 In order to perform to *Modal Analysis* and to obtain first a *Response Model*, the following devices are used:
 - An *acquisition system* (OROS) with 24bits ADCs (figure [[fig:oros]])
@@ -188,38 +201,51 @@ The precise determination of the position of each accelerometer is done using th
 
 The precise position of all the 23 accelerometer with respect to a frame located at the point of interest (located 270mm above the top platform of the hexapod) is shown below. The values are in meter.
 They are contained in the =mat/id31_nanostation.cfg= file.
-
-#+begin_src bash :exports results :eval no-export :post addhdr(*this*)
-  echo " ID X[m] Y[m] Z[m]"
-  cat mat/id31_nanostation.cfg | grep NODES -A 23 | sed '/\s\+[0-9]\+/!d' | sed 's/\(.*\)\s\+0\s\+.\+/\1/' | tac --
+#+begin_src bash :results none
+  cat mat/id31_nanostation.cfg | grep NODES -A 23 | sed '/\s\+[0-9]\+/!d' | sed 's/\(.*\)\s\+0\s\+.\+/\1/' > mat/acc_pos.txt
 #+end_src
 
+We then import that on =matlab=, and sort them.
+#+begin_src matlab
+  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);
+#+end_src
+
+The positions of the sensors relative to the point of interest are shown below (table [[tab:position_accelerometers]]).
+#+begin_src matlab :exports results :results value table replace :post addhdr(*this*)
+  data2orgtable([[1:23]', 1000*acc_pos], {}, {'ID', 'x [mm]', 'y [mm]', 'z [mm]'}, ' %.0f ');
+#+end_src
+
+#+name: tab:position_accelerometers
+#+caption: position of the accelerometers
 #+RESULTS:
-| ID |         X[m] |         Y[m] |         Z[m] |
-|----+--------------+--------------+--------------|
-|  1 | -6.4000e-002 | -6.4000e-002 | -2.7000e-001 |
-|  2 | -6.4000e-002 |  6.4000e-002 | -2.7000e-001 |
-|  3 |  6.4000e-002 |  6.4000e-002 | -2.7000e-001 |
-|  4 |  6.4000e-002 | -6.4000e-002 | -2.7000e-001 |
-|  5 | -3.8500e-001 | -3.0000e-001 | -4.1680e-001 |
-|  6 | -4.2000e-001 |  2.8000e-001 | -4.1680e-001 |
-|  7 |  4.2000e-001 |  2.8000e-001 | -4.1680e-001 |
-|  8 |  3.8000e-001 | -3.0000e-001 | -4.1680e-001 |
-|  9 | -4.7500e-001 | -4.1400e-001 | -4.2730e-001 |
-| 10 | -4.6500e-001 |  4.0700e-001 | -4.2730e-001 |
-| 11 |  4.7500e-001 |  4.2400e-001 | -4.2730e-001 |
-| 12 |  4.7500e-001 | -4.1900e-001 | -4.2730e-001 |
-| 13 | -3.2000e-001 | -4.4600e-001 | -7.8560e-001 |
-| 14 | -4.8000e-001 |  5.3400e-001 | -7.8560e-001 |
-| 15 |  4.5000e-001 |  5.3400e-001 | -7.8560e-001 |
-| 16 |  2.9500e-001 | -4.8100e-001 | -7.8560e-001 |
-| 17 | -7.3000e-001 | -5.2600e-001 | -9.5060e-001 |
-| 18 | -7.3500e-001 |  8.1400e-001 | -9.5060e-001 |
-| 19 |  8.7500e-001 |  7.9900e-001 | -9.5060e-001 |
-| 20 |  8.6490e-001 | -5.0600e-001 | -9.5060e-001 |
-| 21 | -1.5500e-001 | -9.0000e-002 | -5.9400e-001 |
-| 22 |  0.0000e+000 |  1.8000e-001 | -5.9400e-001 |
-| 23 |  1.5500e-001 | -9.0000e-002 | -5.9400e-001 |
+| ID | x [mm] | y [mm] | z [mm] |
+|----+--------+--------+--------|
+|  1 |    -64 |    -64 |   -270 |
+|  2 |    -64 |     64 |   -270 |
+|  3 |     64 |     64 |   -270 |
+|  4 |     64 |    -64 |   -270 |
+|  5 |   -385 |   -300 |   -417 |
+|  6 |   -420 |    280 |   -417 |
+|  7 |    420 |    280 |   -417 |
+|  8 |    380 |   -300 |   -417 |
+|  9 |   -475 |   -414 |   -427 |
+| 10 |   -465 |    407 |   -427 |
+| 11 |    475 |    424 |   -427 |
+| 12 |    475 |   -419 |   -427 |
+| 13 |   -320 |   -446 |   -786 |
+| 14 |   -480 |    534 |   -786 |
+| 15 |    450 |    534 |   -786 |
+| 16 |    295 |   -481 |   -786 |
+| 17 |   -730 |   -526 |   -951 |
+| 18 |   -735 |    814 |   -951 |
+| 19 |    875 |    799 |   -951 |
+| 20 |    865 |   -506 |   -951 |
+| 21 |   -155 |    -90 |   -594 |
+| 22 |      0 |    180 |   -594 |
+| 23 |    155 |    -90 |   -594 |
 
 ** Hammer Impacts
 Only 3 impact points are used.
@@ -635,9 +661,6 @@ We generate such FRF matrix from the measurements using the following script.
         meas_factor = meas_factor*(-1);
       end
 
-      % FRFs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
-      % COHs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
-
       FRFs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
       COHs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
     end
@@ -649,3 +672,221 @@ And we save the obtained FRF matrix and Coherence matrix in a =.mat= file.
 #+begin_src matlab
   save('./mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
 #+end_src
+* Plot showing the coherence of all the measurements
+Now that we have defined a Coherence matrix, we can plot each of its elements to have an idea of the overall coherence and thus, quality of the measurement.
+The result is shown on figure [[fig:all_coherence]].
+
+#+begin_src matlab :exports none
+  n_acc = 23;
+
+  figure;
+  hold on;
+  for i = 1:3*n_acc
+    plot(freqs, squeeze(COHs(i, 1, :)), 'color', [0, 0, 0, 0.2]);
+  end
+  hold off;
+  xlabel('Frequency [Hz]');
+  ylabel('Coherence [\%]');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/all_coherence.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:all_coherence
+#+CAPTION: Plot of the coherence of all the measurements
+[[file:figs/all_coherence.png]]
+
+* Solid Bodies considered for further analysis
+We consider the following solid bodies for further analysis:
+- Bottom Granite
+- Top Granite
+- Translation Stage
+- Tilt Stage
+- Spindle
+- Hexapod
+
+We create a =matlab= structure =solids= that contains the accelerometers ID connected to each solid bodies (as shown on figure [[fig:nass-modal-test]]).
+#+begin_src matlab
+  solids = {};
+  solids.gbot = [17, 18, 19, 20];
+  solids.gtop = [13, 14, 15, 16];
+  solids.ty   = [9, 10, 11, 12];
+  solids.ry   = [5, 6, 7, 8];
+  solids.rz   = [21, 22, 23];
+  solids.hexa = [1, 2, 3, 4];
+
+  solid_names = fields(solids);
+#+end_src
+
+Finally, we save that into a =.mat= file.
+#+begin_src matlab
+  save('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
+#+end_src
+
+* Note about the solid body assumption
+If we measure the motion of a rigid body along a direction $\vec{x}$ using 2 sensors that are co-linear with the same direction $\vec{x}$ ($\vec{p}_2 = \vec{p}_1 + \alpha \vec{x}$), they will measured the same quantity.
+
+This is illustrated on figure [[fig:aligned_accelerometers]].
+
+#+begin_src latex :file aligned_accelerometers.pdf :post pdf2svg(file=*this*, ext="png") :exports results
+  \newcommand\irregularcircle[2]{% radius, irregularity
+    \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
+    +(0:\len pt)
+    \foreach \a in {10,20,...,350}{
+      \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
+      -- +(\a:\len pt)
+    } -- cycle
+  }
+  \begin{tikzpicture}
+    \draw[rounded corners=1mm, fill=blue!30!white] (0, 0) \irregularcircle{3cm}{1mm};
+
+    \node[] (origin) at (4, -1) {$\bullet$};
+    \begin{scope}[shift={(origin)}]
+      \def\axissize{0.8cm}
+      \draw[->] (0, 0) -- ++(\axissize, 0) node[above left]{$x$};
+      \draw[->] (0, 0) -- ++(0, \axissize) node[below right]{$y$};
+      \draw[fill, color=black] (0, 0) circle (0.05*\axissize);
+      \node[draw, circle, inner sep=0pt, minimum size=0.4*\axissize, label=left:$z$] (yaxis) at (0, 0){};
+      \node[below right] at (0, 0){$\{O\}$};
+    \end{scope}
+
+    \coordinate[] (p1) at (-1.5,  1.5);
+    \coordinate[] (p2) at ( 1.5,  1.5);
+
+    \draw[->] (p1)node[]{$\bullet$}node[above]{$p_1$} -- ++(1, 0)node[above]{$v_{x1}$};
+    \draw[->] (p2)node[]{$\bullet$}node[above]{$p_2$} -- ++(1, 0)node[above]{$v_{x2}$};
+
+    \draw[dashed] ($(p1)+(-1, 0)$) -- ($(p2)+(2, 0)$);
+  \end{tikzpicture}
+#+end_src
+
+#+name: fig:aligned_accelerometers
+#+caption: Aligned measurement of the motion of a solid body
+#+RESULTS:
+[[file:figs/aligned_accelerometers.png]]
+
+The motion of the rigid body of figure [[fig:aligned_accelerometers]] is defined by the velocity $\vec{v}$ and rotation $\vec{\Omega}$ with respect to the reference frame $\{O\}$.
+
+The motions at points $1$ and $2$ are:
+\begin{align*}
+  v_1 &= v + \Omega \times p_1 \\
+  v_2 &= v + \Omega \times p_2
+\end{align*}
+
+Taking only the $x$ direction:
+\begin{align*}
+  v_{x1} &= v + \Omega_y p_{z1} - \Omega_z p_{y1} \\
+  v_{x2} &= v + \Omega_y p_{z2} - \Omega_z p_{y2}
+\end{align*}
+
+However, we have $p_{1y} = p_{2y}$ and $p_{1z} = p_{2z}$ because of the co-linearity of the two sensors in the $x$ direction, and thus we obtain
+\begin{equation}
+  v_{x1} = v_{x2}
+\end{equation}
+
+#+begin_important
+  Two sensors that are measuring the motion of a rigid body in the direction of the line linking the two sensors should measure the same quantity.
+#+end_important
+
+We can verify that the rigid body assumption is correct by comparing the measurement of the sensors.
+
+From the table [[tab:position_accelerometers]], we can guess which sensors will give the same results in the X and Y directions.
+
+Comparison of such measurements in the X direction is shown on figure [[fig:compare_acc_x_dir]] and in the Y direction on figure [[fig:compare_acc_y_dir]].
+
+#+begin_src matlab :exports none
+  meas_dir = 1;
+  exc_dir = 1;
+
+  acc_i = [1 , 4 ;
+           2 , 3 ;
+           5 , 8 ;
+           6 , 7 ;
+           9 , 12;
+           10, 11;
+           14, 15;
+           18, 19;
+           21, 23];
+
+  figure;
+  for i = 1:size(acc_i, 1)
+    subaxis(3, 3, i);
+    hold on;
+    plot(freqs, abs(squeeze(FRFs(meas_dir+3*(acc_i(i, 1)-1), exc_dir, :))))
+    plot(freqs, abs(squeeze(FRFs(meas_dir+3*(acc_i(i, 2)-1), exc_dir, :))))
+    hold off;
+    set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+    if i > 6
+      xlabel('Frequency [Hz]');
+    else
+      set(gca, 'XTickLabel',[]);
+    end
+
+    if rem(i, 3) == 1
+      ylabel('Amplitude');
+    end
+    xlim([1, 200]);
+    title(sprintf('Acc %i and %i - X', acc_i(i, 1), acc_i(i, 2)));
+  end
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/compare_acc_x_dir.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:compare_acc_x_dir
+#+CAPTION: Compare accelerometers align in the X direction
+[[file:figs/compare_acc_x_dir.png]]
+
+
+#+begin_src matlab :exports none
+  meas_dir = 2;
+  exc_dir = 1;
+
+  acc_i = [1, 2;
+           5, 6;
+           7, 8;
+           9, 10;
+           11, 12;
+           13, 14;
+           15, 16;
+           17, 18;
+           19, 20];
+
+  figure;
+  for i = 1:size(acc_i, 1)
+    subaxis(3, 3, i);
+    hold on;
+    plot(freqs, abs(squeeze(FRFs(meas_dir+3*(acc_i(i, 1)-1), exc_dir, :))))
+    plot(freqs, abs(squeeze(FRFs(meas_dir+3*(acc_i(i, 2)-1), exc_dir, :))))
+    hold off;
+    set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+    if i > 6
+      xlabel('Frequency [Hz]');
+    else
+      set(gca, 'XTickLabel',[]);
+    end
+
+    if rem(i, 3) == 1
+      ylabel('Amplitude');
+    end
+    xlim([1, 200]);
+    title(sprintf('Acc %i and %i - Y', acc_i(i, 1), acc_i(i, 2)));
+  end
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/compare_acc_y_dir.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:compare_acc_y_dir
+#+CAPTION: Compare accelerometers align in the Y direction
+[[file:figs/compare_acc_y_dir.png]]
+
+#+begin_important
+  From the two figures above, we are more confident about the rigid body assumption in the frequency band of interest.
+#+end_important

modal-analysis/modal_extraction.org

@@ -48,28 +48,10 @@
   <<matlab-init>>
 #+end_src
 
-
-
-* TODO After that, this should be in the modal extraction
-* TODO Part to explain how to choose the modes
-* Obtained Mode Shapes animations
-One all the FRFs are obtained, we can estimate the modal parameters (resonance frequencies, modal shapes and modal damping) within the modal software.
-
-For that, multiple modal extraction techniques can be used (SIMO, MIMO, narrow band, wide band, ...).
-
-Then, it is possible to show the modal shapes with an animation.
-
-Examples are shown on figures [[fig:mode1]] and [[fig:mode6]].
-
-Animations of all the other modes are accessible using the following links: [[file:img/modes/mode1.gif][mode 1]], [[file:img/modes/mode2.gif][mode 2]], [[file:img/modes/mode3.gif][mode 3]], [[file:img/modes/mode4.gif][mode 4]], [[file:img/modes/mode5.gif][mode 5]], [[file:img/modes/mode6.gif][mode 6]], [[file:img/modes/mode7.gif][mode 7]], [[file:img/modes/mode8.gif][mode 8]], [[file:img/modes/mode9.gif][mode 9]], [[file:img/modes/mode10.gif][mode 10]], [[file:img/modes/mode11.gif][mode 11]], [[file:img/modes/mode12.gif][mode 12]], [[file:img/modes/mode13.gif][mode 13]], [[file:img/modes/mode14.gif][mode 14]], [[file:img/modes/mode15.gif][mode 15]], [[file:img/modes/mode16.gif][mode 16]], [[file:img/modes/mode17.gif][mode 17]], [[file:img/modes/mode18.gif][mode 18]], [[file:img/modes/mode19.gif][mode 19]], [[file:img/modes/mode20.gif][mode 20]], [[file:img/modes/mode21.gif][mode 21]].
-
-#+name: fig:mode1
-#+caption: Mode 1
-[[file:img/modes/mode1.gif]]
-
-#+name: fig:mode6
-#+caption: Mode 6
-[[file:img/modes/mode6.gif]]
+* TODO Part to explain how to choose the modes frequencies
+- bro-band method used
+- Stabilization Chart
+- 21 modes
 
 * Obtained Modal Parameters
 From the modal analysis software, we can export the obtained modal parameters:
@@ -212,6 +194,25 @@ The obtained mode frequencies and damping are shown below.
 |          150.1 |         2.2 |
 |          164.7 |         1.4 |
 
+* Obtained Mode Shapes animations
+One all the FRFs are obtained, we can estimate the modal parameters (resonance frequencies, modal shapes and modal damping) within the modal software.
+
+For that, multiple modal extraction techniques can be used (SIMO, MIMO, narrow band, wide band, ...).
+
+Then, it is possible to show the modal shapes with an animation.
+
+Examples are shown on figures [[fig:mode1]] and [[fig:mode6]].
+
+Animations of all the other modes are accessible using the following links: [[file:img/modes/mode1.gif][mode 1]], [[file:img/modes/mode2.gif][mode 2]], [[file:img/modes/mode3.gif][mode 3]], [[file:img/modes/mode4.gif][mode 4]], [[file:img/modes/mode5.gif][mode 5]], [[file:img/modes/mode6.gif][mode 6]], [[file:img/modes/mode7.gif][mode 7]], [[file:img/modes/mode8.gif][mode 8]], [[file:img/modes/mode9.gif][mode 9]], [[file:img/modes/mode10.gif][mode 10]], [[file:img/modes/mode11.gif][mode 11]], [[file:img/modes/mode12.gif][mode 12]], [[file:img/modes/mode13.gif][mode 13]], [[file:img/modes/mode14.gif][mode 14]], [[file:img/modes/mode15.gif][mode 15]], [[file:img/modes/mode16.gif][mode 16]], [[file:img/modes/mode17.gif][mode 17]], [[file:img/modes/mode18.gif][mode 18]], [[file:img/modes/mode19.gif][mode 19]], [[file:img/modes/mode20.gif][mode 20]], [[file:img/modes/mode21.gif][mode 21]].
+
+#+name: fig:mode1
+#+caption: Mode 1
+[[file:img/modes/mode1.gif]]
+
+#+name: fig:mode6
+#+caption: Mode 6
+[[file:img/modes/mode6.gif]]
+
 * Compute the Modal Model
 ** Matlab Init                                              :noexport:ignore:
 #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@@ -1201,7 +1202,6 @@ This could be due to the 4 Airloc Levelers that are used for the granite (figure
 
 They are probably *not well leveled*, so the granite is supported only by two Airloc.
 
-
 * Setup
 #+name: fig:nass-modal-test
 #+caption: Position and orientation of the accelerometer used
@@ -1356,7 +1356,7 @@ We create a structure =solids= that contains the accelerometer number of each so
   solid_names = fields(solids);
 #+end_src
 
-* From local coordinates to global coordinates
+* From local coordinates to global coordinates for the mode shapes
 #+begin_src latex :file local_to_global_coordinates.pdf :post pdf2svg(file=*this*, ext="png") :exports results
   \newcommand\irregularcircle[2]{% radius, irregularity
     \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
@@ -1791,21 +1791,7 @@ FRF matrix $n \times p$:
   xlim([1, 200]);
 #+end_src
 
-
-#+begin_src matlab
-  figure;
-  hold on;
-  for i = 1:3*n_acc
-    plot(freqs, squeeze(COHs(i, 1, :)), 'color', [0, 0, 0, 0.2]);
-  end
-  hold off;
-  xlabel('Frequency [Hz]');
-  ylabel('Coherence [\%]');
-#+end_src
-
-
-Composite Response Function.
-
+* Composite Response Function
 We here sum the norm instead of the complex numbers.
 
 #+begin_src matlab
@@ -1833,6 +1819,30 @@ We here sum the norm instead of the complex numbers.
   xlim([1, 200]);
 #+end_src
 
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/composite_response_function.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:composite_response_function
+#+CAPTION: Composite Response Function
+[[file:figs/composite_response_function.png]]
+
+
+* TODO Singular Value Decomposition - Modal Indication Function
+Show the same plot as in the modal software.
+This helps to identify double modes.
+
+From the documentation of the modal software:
+#+begin_quote
+The MIF consist of the singular values of the Frequency response function matrix.
+The number of MIFs equals the number of excitations.
+By the powerful singular value decomposition, the real signal space is separated from the noise space.
+Therefore, the MIFs exhibit the modes effectively.
+A peak in the MIFs plot usually indicate the existence of a structural mode, and two peaks at the same frequency point means the existence of two repeated modes.
+Moreover, the magnitude of the MIFs implies the strength of the a mode.
+#+end_quote
+
 * From local coordinates to global coordinates with the FRFs
 #+begin_src matlab
   % Number of Solids * DOF for each solid / Number of excitation / frequency points