Add many data analysis. Clean modal_extraction file

modal-analysis/frf_processing.org

@@ -125,10 +125,10 @@ The goal here is to link these $4 \times 3 = 12$ measurements to the 6 DOFs of t
     \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$};
+    \draw[->] (p1)node[]{$\bullet$}node[above]{$p_1$} -- ++( 1  , 0.5)node[right]{$\delta p_1$};
+    \draw[->] (p2)node[]{$\bullet$}node[above]{$p_2$} -- ++(-0.5, 1  )node[right]{$\delta p_2$};
+    \draw[->] (p3)node[]{$\bullet$}node[above]{$p_3$} -- ++( 1  , 0.5)node[right]{$\delta p_3$};
+    \draw[->] (p4)node[]{$\bullet$}node[above]{$p_4$} -- ++( 0.5, 1  )node[right]{$\delta p_4$};
   \end{tikzpicture}
 #+end_src
 
@@ -137,23 +137,23 @@ The goal here is to link these $4 \times 3 = 12$ measurements to the 6 DOFs of t
 #+RESULTS:
 [[file:figs/local_to_global_coordinates.png]]
 
+We consider the motion of the rigid body defined by its displacement $\delta p$ and rotation $\vec{\Omega}$ with respect to the reference frame $\{O\}$.
+
 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
+  \delta p_1 &= \delta p + \delta\Omega p_1\\
+  \delta p_2 &= \delta p + \delta\Omega p_2\\
+  \delta p_3 &= \delta p + \delta\Omega p_3\\
+  \delta p_4 &= \delta p + \delta\Omega 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
+  \delta\Omega = \begin{bmatrix}
+   0              & -\delta\Omega_z &  \delta\Omega_y \\
+   \delta\Omega_z &  0              & -\delta\Omega_x \\
+  -\delta\Omega_y &  \delta\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:
 
@@ -167,16 +167,16 @@ We can rearrange the equations in a matrix form:
   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
+  \delta p_x \\ \delta p_y \\ \delta p_z \\ \hline \delta\Omega_x \\ \delta\Omega_y \\ \delta\Omega_z
 \end{bmatrix} = \begin{bmatrix}
-  v_{1x} \\ v_{1y} \\ v_{1z} \\\hline \vdots \\\hline v_{4x} \\ v_{4y} \\ v_{4z}
+  \delta p_{1x} \\ \delta p_{1y} \\ \delta p_{1z} \\\hline \vdots \\\hline \delta p_{4x} \\ \delta p_{4y} \\ \delta p_{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
+  \delta p_x \\ \delta p_y \\ \delta p_z \\ \hline \delta\Omega_x \\ \delta\Omega_y \\ \delta\Omega_z
 \end{bmatrix} =
 \left[\begin{array}{ccc|ccc}
   1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\
@@ -187,18 +187,18 @@ and then we obtain the velocity and rotation of the solid in the wanted frame $\
   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}
+  \delta p_{1x} \\ \delta p_{1y} \\ \delta p_{1z} \\\hline \vdots \\\hline \delta p_{4x} \\ \delta p_{4y} \\ \delta p_{4z}
+\end{bmatrix} \label{eq:determine_global_disp}
 \end{equation}
 
-This inversion is equivalent to resolving a mean square problem.
+#+begin_important
+  Using equation eqref:eq:determine_global_disp, we can determine the motion of the solid body expressed in a chosen frame $\{O\}$ using the accelerometers attached to it.
+  The inversion is equivalent to resolving a mean square problem.
+#+end_important
 
 * 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)
@@ -212,9 +212,15 @@ The possibles choices are:
 | 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.
+The choice of the frame depends of what we want to do with the data.
 
-* From accelerometer DOFs to solid body DOFs - Matlab Implementation
+One of the goals 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.
+
+In the next sections, we will express the FRF matrix in the different frames.
+
+* FRF expressed in a common frame
+** 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$
@@ -223,7 +229,7 @@ First, we initialize a new FRF matrix =FRFs_O= which is an $n \times p \times q$
 #+begin_important
 For each frequency point $\omega_i$, the FRF matrix =FRFs_O= is a $n\times p$ matrix:
 \begin{equation}
-  \text{FRF}_O(\omega_i) = \begin{bmatrix}
+  \text{FRF}_\text{O}(\omega_i) = \begin{bmatrix}
   \frac{D_{1,T_x}}{F_x}(\omega_i) & \frac{D_{1,T_x}}{F_y}(\omega_i) & \frac{D_{1,T_x}}{F_z}(\omega_i) \\
   \frac{D_{1,T_y}}{F_x}(\omega_i) & \frac{D_{1,T_y}}{F_y}(\omega_i) & \frac{D_{1,T_y}}{F_z}(\omega_i) \\
   \frac{D_{1,T_z}}{F_x}(\omega_i) & \frac{D_{1,T_z}}{F_y}(\omega_i) & \frac{D_{1,T_z}}{F_z}(\omega_i) \\
@@ -235,15 +241,14 @@ For each frequency point $\omega_i$, the FRF matrix =FRFs_O= is a $n\times p$ ma
   \frac{D_{6,R_z}}{F_x}(\omega_i) & \frac{D_{6,R_z}}{F_y}(\omega_i) & \frac{D_{6,R_z}}{F_z}(\omega_i)
 \end{bmatrix}
 \end{equation}
-where 1, 2, ..., 6 corresponds to the 6 solid bodies.
+where $D_i$ corresponds to the solid body number i.
 #+end_important
 
+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
   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});
 
@@ -263,7 +268,7 @@ Then, as we know the positions of the accelerometers on each solid body, and we
   end
 #+end_src
 
-* Analysis of some FRF in the global coordinates
+** 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]]).
@@ -285,7 +290,7 @@ We can also compare all the DOFs of one solid body (figure [[fig:frf_one_body_al
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
-  ylabel('Amplitude');
+  ylabel('Amplitude [$\frac{m/s^2}{N}$]');
   legend('Location', 'northwest');
   title(sprintf('FRF between %s and %s', exc_names{exc_dir}, DOFs{dir_i}));
 
@@ -328,7 +333,7 @@ We can also compare all the DOFs of one solid body (figure [[fig:frf_one_body_al
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   set(gca, 'XTickLabel',[]);
-  ylabel('Amplitude');
+  ylabel('Amplitude [$\frac{m/s^2}{N}$]');
   legend('Location', 'northwest');
   title(sprintf('Motion of %s due to %s', solid_names{solid_i}, exc_names{exc_dir}));
 
@@ -355,15 +360,15 @@ We can also compare all the DOFs of one solid body (figure [[fig:frf_one_body_al
 #+CAPTION: FRFs of one solid body in all its DOFs
 [[file:figs/frf_one_body_all_directions.png]]
 
-* Comparison of the relative motion of solid bodies
+** Comparison of the relative motion of solid bodies
 Now that the motion of all the solid bodies are expressed in the same frame, we should be able to *compare them*.
 This can be used to determine what joints direction between two solid bodies is stiff enough that we can fix this DoF.
 This could help reduce the order of the model and simplify the extraction of the model parameters from the measurements.
 
 We decide to plot the "normalized relative motion" between solid bodies $i$ and $j$:
-\[ 0 < \Delta_{ij, x} = \frac{\left| D_{i,x} - D_{j,x} \right|}{|D_{i,x}| + |D_{j,x}|} < 1 \]
+\[ 0 < \Delta_{ij, x} = \frac{\left| |D_{i,x}| - |D_{j,x}| \right|}{|D_{i,x}| + |D_{j,x}|} < 1 \]
 
-Then, if $\Delta_{ij,x} \ll 0$ in the frequency band of interest, we have that $D_{ix} \approx D_{jx}$ and we can neglect that DOF between the two solid bodies $i$ and $j$.
+Then, if $\Delta_{ij,x} \ll 1$ in the frequency band of interest, we have that $D_{ix} \approx D_{jx}$ and we can neglect that DOF between the two solid bodies $i$ and $j$.
 
 This normalized relative motion is shown on figure [[fig:relative_motion_comparison]] for all the directions and for all the adjacent pair of solid bodies.
 
@@ -379,7 +384,7 @@ This normalized relative motion is shown on figure [[fig:relative_motion_compari
     subplot(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, :))));
+      H = (abs(squeeze(FRFs_O((i-1)*6+dir_i, exc_dir, :)))-abs(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;
@@ -420,18 +425,29 @@ This normalized relative motion is shown on figure [[fig:relative_motion_compari
   The relative motion of two solid bodies may be negligible when exciting the structure at on point and but at another point.
 #+end_warning
 
-* 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.
+** 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 with respect to the frame $\{O\}$.
 
-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.
+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, supposing that the stage is a solid body.
 
 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
 
+From the translation $\delta p$ and rotation $\delta \Omega$ of a solid body and the positions $p_i$ of the accelerometers attached to it, we can compute the response that would have been measured by the accelerometers using the following formula:
+\begin{align*}
+  \delta p_1 &= \delta p + \delta\Omega p_1\\
+  \delta p_2 &= \delta p + \delta\Omega p_2\\
+  \delta p_3 &= \delta p + \delta\Omega p_3\\
+  \delta p_4 &= \delta p + \delta\Omega p_4
+\end{align*}
+
+Thus, we can obtain the FRF matrix =FRFs_A= that gives the responses of the accelerometers to the forces applied by the hammer.
+
+It is implemented in matlab as follow:
 #+begin_src matlab
-  FRF_recovered = zeros(size(FRFs));
+  FRFs_A = zeros(size(FRFs));
 
   % For each excitation direction
   for exc_dir = 1:3
@@ -443,16 +459,16 @@ This will help us to determine if:
       % 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, :)';
+        pos = acc_pos(acc_i, :).';
         posX = [0 pos(3) -pos(2); -pos(3) 0 pos(1) ; pos(2) -pos(1) 0];
 
-        FRF_recovered(3*(acc_i-1)+1:3*(acc_i-1)+3, exc_dir, :) = v0 + posX*W0;
+        FRFs_A(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.
+We then compare the original FRF measured for each accelerometer =FRFs= with the "recovered" FRF =FRFs_A= 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.
@@ -482,7 +498,7 @@ The FRF are matching well until 100Hz.
     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');
+      plot(freqs, abs(squeeze(FRFs_A(3*(acc_i-1)+dir_i, exc_dir, :))), '--', 'HandleVisibility', 'off');
     end
     hold off;
 
@@ -494,10 +510,10 @@ The FRF are matching well until 100Hz.
     end
 
     if rem(i, 2) == 1
-      ylabel('Amplitude');
+      ylabel('Amplitude [$\frac{m/s^2}{N}$]');
     end
 
-    xlim([1, 200]);
+    xlim([1, 200]); ylim([1e-6, 1e-1]);
     title(sprintf('Accelerometer %i', accs_i(i)));
     legend('location', 'northwest');
   end
@@ -534,7 +550,7 @@ The FRF are matching well until 100Hz.
     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');
+      plot(freqs, abs(squeeze(FRFs_A(3*(acc_i-1)+dir_i, exc_dir, :))), '--', 'HandleVisibility', 'off');
     end
     hold off;
 
@@ -546,10 +562,10 @@ The FRF are matching well until 100Hz.
     end
 
     if rem(i, 2) == 1
-      ylabel('Amplitude');
+      ylabel('Amplitude [$\frac{m/s^2}{N}$]');
     end
 
-    xlim([1, 200]);
+    xlim([1, 200]); ylim([1e-6, 1e-1]);
     title(sprintf('Accelerometer %i', accs_i(i)));
     legend('location', 'northwest');
   end
@@ -572,7 +588,187 @@ The FRF are matching well until 100Hz.
   This valid the fact that a multi-body model can be used to represent the dynamics of the micro-station.
 #+end_important
 
-* Saving of the FRF expressed in the global coordinates
+** Saving of the FRF expressed in the global coordinates
 #+begin_src matlab
   save('mat/frf_o.mat', 'FRFs_O');
 #+end_src
+
+* FRF expressed in a frame centered at the CoM of each solid body
+** Center of Mass of each solid body
+From solidworks, we can export the position of the center of mass of each solid body considered.
+
+For that, we have to suppose density of the materials.
+We then process the exported file.
+
+#+begin_src bash :results none
+  sed -e '1,/^Center of mass/ d' mat/1_granite_bot.txt | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >  mat/model_solidworks_com.txt
+  sed -e '1,/^Center of mass/ d' mat/2_granite_top.txt | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >> mat/model_solidworks_com.txt
+  sed -e '1,/^Center of mass/ d' mat/3_ty.txt          | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >> mat/model_solidworks_com.txt
+  sed -e '1,/^Center of mass/ d' mat/4_ry.txt          | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >> mat/model_solidworks_com.txt
+  sed -e '1,/^Center of mass/ d' mat/5_rz.txt          | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >> mat/model_solidworks_com.txt
+  sed -e '1,/^Center of mass/ d' mat/6_hexa.txt        | sed 3q | sed 's/^\s*\t*[XYZ] = \([-0-9.]*\)\r/\1/' >> mat/model_solidworks_com.txt
+#+end_src
+
+And we import the Center of mass on Matlab.
+#+begin_src matlab
+  model_com = table2array(readtable('mat/model_solidworks_com.txt', 'ReadVariableNames', false));
+#+end_src
+
+#+begin_src matlab
+  model_com = reshape(model_com , [3, 6]);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable(1000*model_com, {'X [mm]', 'Y [mm]', 'Z [mm]'}, {'granite bot', 'granite top', 'ty', 'ry', 'rz', 'hexa'}, ' %.0f ');
+#+end_src
+
+#+RESULTS:
+|        | granite bot | granite top |   ty |   ry |   rz | hexa |
+|--------+-------------+-------------+------+------+------+------|
+| X [mm] |          45 |          52 |    0 |    0 |    0 |   -4 |
+| Y [mm] |         144 |         258 |   14 |   -5 |    0 |    6 |
+| Z [mm] |       -1251 |        -778 | -600 | -628 | -580 | -319 |
+
+** From accelerometer DOFs to solid body DOFs - Expressed at the CoM
+First, we initialize a new FRF matrix =FRFs_CoM= 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_important
+For each frequency point $\omega_i$, the FRF matrix =FRFs_CoM= is a $n\times p$ matrix:
+\begin{equation}
+  \text{FRF}_\text{CoM}(\omega_i) = \begin{bmatrix}
+  \frac{D_{1,T_x}}{F_x}(\omega_i) & \frac{D_{1,T_x}}{F_y}(\omega_i) & \frac{D_{1,T_x}}{F_z}(\omega_i) \\
+  \frac{D_{1,T_y}}{F_x}(\omega_i) & \frac{D_{1,T_y}}{F_y}(\omega_i) & \frac{D_{1,T_y}}{F_z}(\omega_i) \\
+  \frac{D_{1,T_z}}{F_x}(\omega_i) & \frac{D_{1,T_z}}{F_y}(\omega_i) & \frac{D_{1,T_z}}{F_z}(\omega_i) \\
+  \frac{D_{1,R_x}}{F_x}(\omega_i) & \frac{D_{1,R_x}}{F_y}(\omega_i) & \frac{D_{1,R_x}}{F_z}(\omega_i) \\
+  \frac{D_{1,R_y}}{F_x}(\omega_i) & \frac{D_{1,R_y}}{F_y}(\omega_i) & \frac{D_{1,R_y}}{F_z}(\omega_i) \\
+  \frac{D_{1,R_z}}{F_x}(\omega_i) & \frac{D_{1,R_z}}{F_y}(\omega_i) & \frac{D_{1,R_z}}{F_z}(\omega_i) \\
+  \frac{D_{2,T_x}}{F_x}(\omega_i) & \frac{D_{2,T_x}}{F_y}(\omega_i) & \frac{D_{2,T_x}}{F_z}(\omega_i) \\
+  \vdots              & \vdots              & \vdots              \\
+  \frac{D_{6,R_z}}{F_x}(\omega_i) & \frac{D_{6,R_z}}{F_y}(\omega_i) & \frac{D_{6,R_z}}{F_z}(\omega_i)
+\end{bmatrix}
+\end{equation}
+where 1, 2, ..., 6 corresponds to the 6 solid bodies.
+#+end_important
+
+#+begin_src matlab
+  FRFs_CoM = 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);
+
+      acc_pos_com = acc_pos(acc_i, :).' - model_com(:, solid_i);
+
+      A(3*(i-1)+1:3*i, 1:3) = eye(3);
+      A(3*(i-1)+1:3*i, 4:6) = [ 0               acc_pos_com(3) -acc_pos_com(2) ;
+                               -acc_pos_com(3)  0               acc_pos_com(1) ;
+                                acc_pos_com(2) -acc_pos_com(1)  0];
+    end
+
+    for exc_dir = 1:3
+      FRFs_CoM((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
+#+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 = subplot(2, 1, 1);
+  hold on;
+  for solid_i = solids_i
+    plot(freqs, abs(squeeze(FRFs_CoM((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 [$\frac{m/s^2}{N}$]');
+  legend('Location', 'northwest');
+  title(sprintf('FRF between %s and %s', exc_names{exc_dir}, DOFs{dir_i}));
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for solid_i = solids_i
+    plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_CoM((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_com_all_bodies_one_direction.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:frf_com_all_bodies_one_direction
+#+CAPTION: FRFs of all the 6 solid bodies in one direction
+[[file:figs/frf_com_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 = 6;
+  dirs_i = 1:6;
+  exc_dir = 3;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for dir_i = dirs_i
+    plot(freqs, abs(squeeze(FRFs_CoM((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 [$\frac{m/s^2}{N}$]');
+  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_CoM((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_com_one_body_all_directions.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:frf_com_one_body_all_directions
+#+CAPTION: FRFs of one solid body in all its DOFs (expressed with a frame centered with its center of mass)
+[[file:figs/frf_com_one_body_all_directions.png]]
+
+** Save the FRF
+#+begin_src matlab
+  save('mat/frf_com.mat', 'FRFs_CoM');
+#+end_src

modal-analysis/mat/1_granite_bot.txt

@@ -0,0 +1,35 @@
+Mass properties of Granite-bot
+     Configuration: Default
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+Density = 2700.000 kilograms per cubic meter
+
+Mass = 4860.000 kilograms
+
+Volume = 1.800 cubic meters
+
+Surface area = 10.200 square meters
+
+Center of mass: ( meters )
+	X = 0.045
+	Y = 0.144
+	Z = -1.251
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 1.000,  0.000,  0.000)   	Px = 1057.050
+	 Iy = ( 0.000,  1.000,  0.000)   	Py = 1765.800
+	 Iz = ( 0.000,  0.000,  1.000)   	Pz = 2531.250
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 1057.050	Lxy = 0.000	Lxz = 0.000
+	Lyx = 0.000	Lyy = 1765.800	Lyz = 0.000
+	Lzx = 0.000	Lzy = 0.000	Lzz = 2531.250
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 8758.303	Ixy = 31.499	Ixz = -273.495
+	Iyx = 31.499	Iyy = 9376.075	Iyz = -875.367
+	Izx = -273.495	Izy = -875.367	Izz = 2641.910

modal-analysis/mat/2_granite_top.txt

@@ -0,0 +1,34 @@
+Mass properties of selected components
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+* Includes the mass properties of one or more hidden components/bodies.
+
+Mass = 2152.191 kilograms
+
+Volume = 0.749 cubic meters
+
+Surface area = 14.433 square meters
+
+Center of mass: ( meters )
+	X = 0.052
+	Y = 0.258
+	Z = -0.778
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 1.000, -0.011,  0.000)   	Px = 449.573
+	 Iy = ( 0.011,  1.000,  0.016)   	Py = 1037.927
+	 Iz = ( 0.000, -0.016,  1.000)   	Pz = 1455.367
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 449.640	Lxy = -6.266	Lxz = 0.318
+	Lyx = -6.266	Lyy = 1037.969	Lyz = 6.718
+	Lzx = 0.318	Lzy = 6.718	Lzz = 1455.258
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 1894.930	Ixy = 22.465	Ixz = -86.122
+	Iyx = 22.465	Iyy = 2345.208	Iyz = -425.885
+	Izx = -86.122	Izy = -425.885	Izz = 1604.791

modal-analysis/mat/3_ty.txt

@@ -0,0 +1,34 @@
+Mass properties of selected components
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+* Includes the mass properties of one or more hidden components/bodies.
+
+Mass = 308.383 kilograms
+
+Volume = 0.044 cubic meters
+
+Surface area = 8.343 square meters
+
+Center of mass: ( meters )
+	X = 0.000
+	Y = 0.014
+	Z = -0.600
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 1.000, -0.009,  0.000)   	Px = 40.953
+	 Iy = ( 0.009,  1.000, -0.003)   	Py = 49.447
+	 Iz = ( 0.000,  0.003,  1.000)   	Pz = 84.729
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 40.953	Lxy = -0.074	Lxz = -0.002
+	Lyx = -0.074	Lyy = 49.447	Lyz = -0.105
+	Lzx = -0.002	Lzy = -0.105	Lzz = 84.729
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 151.854	Ixy = -0.073	Ixz = -0.063
+	Iyx = -0.073	Iyy = 160.290	Iyz = -2.637
+	Izx = -0.063	Izy = -2.637	Izz = 84.787

modal-analysis/mat/4_ry.txt

@@ -0,0 +1,34 @@
+Mass properties of selected components
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+* Includes the mass properties of one or more hidden components/bodies.
+
+Mass = 81.987 kilograms
+
+Volume = 0.082 cubic meters
+
+Surface area = 7.022 square meters
+
+Center of mass: ( meters )
+	X = 0.000
+	Y = -0.005
+	Z = -0.628
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 1.000,  0.000,  0.000)   	Px = 5.203
+	 Iy = ( 0.000,  1.000, -0.013)   	Py = 5.575
+	 Iz = ( 0.000,  0.013,  1.000)   	Pz = 7.526
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 5.203	Lxy = 0.000	Lxz = 0.000
+	Lyx = 0.000	Lyy = 5.576	Lyz = -0.025
+	Lzx = 0.000	Lzy = -0.025	Lzz = 7.525
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 37.567	Ixy = 0.000	Ixz = 0.002
+	Iyx = 0.000	Iyy = 37.937	Iyz = 0.241
+	Izx = 0.002	Izy = 0.241	Izz = 7.528

modal-analysis/mat/5_rz.txt

@@ -0,0 +1,34 @@
+Mass properties of selected components
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+* Includes the mass properties of one or more hidden components/bodies.
+
+Mass = 54.689 kilograms
+
+Volume = 0.023 cubic meters
+
+Surface area = 1.873 square meters
+
+Center of mass: ( meters )
+	X = 0.000
+	Y = 0.000
+	Z = -0.580
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 0.707, -0.707,  0.000)   	Px = 1.292
+	 Iy = ( 0.707,  0.707,  0.000)   	Py = 1.292
+	 Iz = ( 0.000,  0.000,  1.000)   	Pz = 1.831
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 1.292	Lxy = 0.000	Lxz = 0.000
+	Lyx = 0.000	Lyy = 1.292	Lyz = 0.000
+	Lzx = 0.000	Lzy = 0.000	Lzz = 1.831
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 19.673	Ixy = 0.000	Ixz = 0.000
+	Iyx = 0.000	Iyy = 19.673	Iyz = 0.000
+	Izx = 0.000	Izy = 0.000	Izz = 1.831

modal-analysis/mat/6_hexa.txt

@@ -0,0 +1,34 @@
+Mass properties of selected components
+     Coordinate system: Coordinate System1
+
+The center of mass and the moments of inertia are output in the coordinate system of Assemblage
+* Includes the mass properties of one or more hidden components/bodies.
+
+Mass = 19.656 kilograms
+
+Volume = 0.005 cubic meters
+
+Surface area = 1.488 square meters
+
+Center of mass: ( meters )
+	X = -0.004
+	Y = 0.006
+	Z = -0.319
+
+Principal axes of inertia and principal moments of inertia: ( kilograms * square meters )
+Taken at the center of mass.
+	 Ix = ( 0.806,  0.593, -0.001)   	Px = 0.250
+	 Iy = (-0.592,  0.805, -0.030)   	Py = 0.361
+	 Iz = (-0.017,  0.024,  1.000)   	Pz = 0.526
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the center of mass and aligned with the output coordinate system.
+	Lxx = 0.289	Lxy = 0.053	Lxz = 0.003
+	Lyx = 0.053	Lyy = 0.322	Lyz = -0.004
+	Lzx = 0.003	Lzy = -0.004	Lzz = 0.526
+
+Moments of inertia: ( kilograms * square meters )
+Taken at the output coordinate system.
+	Ixx = 2.295	Ixy = 0.052	Ixz = 0.031
+	Iyx = 0.052	Iyy = 2.327	Iyz = -0.040
+	Izx = 0.031	Izy = -0.040	Izz = 0.527

modal-analysis/mat/model_solidworks_com.txt

@@ -0,0 +1,18 @@
+0.045
+0.144
+-1.251
+0.052
+0.258
+-0.778
+0.000
+0.014
+-0.600
+0.000
+-0.005
+-0.628
+0.000
+0.000
+-0.580
+-0.004
+0.006
+-0.319

modal-analysis/mat/model_solidworks_inertia.txt

@@ -0,0 +1,13 @@
+8758.303	31.499	-273.495
+31.499	9376.075	-875.367
+-273.495	-875.367	2641.9101894.930	22.465	-86.122
+22.465	2345.208	-425.885
+-86.122	-425.885	1604.791151.854	-0.073	-0.063
+-0.073	160.290	-2.637
+-0.063	-2.637	84.78737.567	0.000	0.002
+0.000	37.937	0.241
+0.002	0.241	7.52819.673	0.000	0.000
+0.000	19.673	0.000
+0.000	0.000	1.8312.295	0.052	0.031
+0.052	2.327	-0.040
+0.031	-0.040	0.527

modal-analysis/mat/model_solidworks_mass.txt

@@ -0,0 +1,6 @@
+4860.000
+2152.191
+308.383
+81.987
+54.689
+19.656

modal-analysis/mathematical_model.org

@@ -61,6 +61,48 @@ For instance, the translation stage may not vibrate in the Z direction for all t
 
 The modal identification done here will thus permit us to determine *which DOF can be neglected*.
 
+* TODO Extract Physical Matrices
+cite:wang11_extrac_real_modes_physic_matric
+
+Let's recall that:
+\[ \Lambda = \begin{bmatrix}
+  s_1 &        & 0 \\
+      & \ddots &   \\
+  0   &        & s_N
+\end{bmatrix}_{N \times N}; \quad \Psi = \begin{bmatrix}
+  & & \\
+  \{\psi_1\} & \dots & \{\psi_N\} \\
+  & &
+\end{bmatrix}_{M \times N} ; \quad A = \begin{bmatrix}
+  a_1 &        & 0   \\
+      & \ddots &     \\
+  0   &        & a_N
+\end{bmatrix}_{N \times N}; \]
+
+
+\begin{align}
+  M &= \frac{1}{2} \left[ \text{Re}(\Psi A^{-1} \Lambda \Psi^T ) \right]^{-1} \\
+  C &= -2 M \text{Re}(\Psi A^{-1} \Lambda^2 A^{-1} \Psi^T ) M \\
+  K &= -\frac{1}{2} \left[ \text{Re}(\Psi \Lambda^{-1} A^{-1} \Psi^T) \right]^{-1}
+\end{align}
+
+#+begin_src matlab
+  psi = eigen_vec_CoM;
+  a = modal_a_M;
+  lambda = eigen_val_M;
+
+  M = 0.5*inv(real(psi*inv(a)*lambda*psi'));
+  C = -2*M*real(psi*inv(a)*lambda^2*inv(a)*psi')*M;
+  K = -0.5*inv(real(psi*inv(lambda)*inv(a)*psi'));
+#+end_src
+
+From cite:ewins00_modal
+
+\begin{align}
+  [M] &= [\Phi]^{-T} [I] [\Phi]^{-1} \\
+  [K] &= [\Phi]^{-T} [\lambda_r^2] [\Phi]^{-1}
+\end{align}
+
 * Some notes about constraining the number of degrees of freedom
 We want to have the two eigen matrices.
 

modal-analysis/matlab/modal_extraction.m

@@ -0,0 +1,768 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+% Load Data
+
+load('mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
+load('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
+load('mat/frf_o.mat', 'FRFs_O');
+
+% Singular Value Decomposition - Modal Indication Function
+% The Mode Indicator Functions are usually used on $n\times p$ FRF matrix where $n$ is a relatively large number of measurement DOFs and $p$ is the number of excitation DOFs, typically 3 or 4.
+
+% In these methods, the frequency dependent FRF matrix is subjected to a singular value decomposition analysis which thus yields a small number (3 or 4) of singular values, these also being frequency dependent.
+
+% These methods are used to *determine the number of modes* present in a given frequency range, to *identify repeated natural frequencies* and to pre process the FRF data prior to modal analysis.
+
+% 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
+
+% #+begin_important
+% The *Complex Mode Indicator Function* is defined simply by the SVD of the FRF (sub) matrix:
+% \begin{align*}
+%   [H(\omega)]_{n\times p} &= [U(\omega)]_{n\times n} [\Sigma(\omega)]_{n\times p} [V(\omega)]_{p\times p}^H\\
+%   [CMIF(\omega)]_{p\times p} &= [\Sigma(\omega)]_{p\times n}^T [\Sigma(\omega)]_{n\times p}
+% \end{align*}
+% #+end_important
+
+% We compute the Complex Mode Indicator Function. The result is shown on figure [[fig:cmif]].
+% The exact same curve is obtained when computed using the OROS software.
+
+
+MIF = zeros(size(FRFs, 2), size(FRFs, 2), size(FRFs, 3));
+
+for i = 1:length(freqs)
+  [~,S,~] = svd(FRFs(:, :, i));
+  MIF(:, :, i) = S'*S;
+end
+
+figure;
+hold on;
+for i = 1:size(MIF, 1)
+  plot(freqs, squeeze(MIF(i, i, :)), 'DisplayName', sprintf('MDIF - %i', i));
+end
+hold off;
+set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('CMIF Amplitude');
+xlim([1, 200]);
+legend('location', 'southeast');
+
+
+
+% #+NAME: fig:cmif
+% #+CAPTION: Complex Mode Indicator Function
+% [[file:figs/cmif.png]]
+
+% We can also compute the CMIF using the FRF matrix expressed in the same global frame.
+% We compare the two CMIF on figure [[fig:cmif_compare]].
+
+% They do not indicate the same resonance frequencies, especially around 110Hz.
+
+
+MIF_O = zeros(size(FRFs_O, 2), size(FRFs_O, 2), size(FRFs_O, 3));
+for i = 1:length(freqs)
+  [~,S,~] = svd(FRFs_O(:, :, i));
+  MIF_O(:, :, i) = S'*S;
+end
+
+figure;
+hold on;
+for i = 1:size(MIF, 1)
+  set(gca,'ColorOrderIndex',i)
+  plot(freqs, squeeze(MIF(i, i, :)), '-');
+  set(gca,'ColorOrderIndex',i)
+  plot(freqs, squeeze(MIF_O(i, i, :)), '--');
+end
+hold off;
+set(gca, 'Yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('CMIF Amplitude');
+xlim([1, 200]);
+
+% Composite Response Function
+% An alternative is the Composite Response Function $HH(\omega)$ defined as the sum of all the measured FRF:
+% \begin{equation}
+%   HH(\omega) = \sum_j\sum_kH_{jk}(\omega)
+% \end{equation}
+
+% Instead, we choose here to use the sum of the norms of the measured FRFs:
+% \begin{equation}
+%   HH(\omega) = \sum_j\sum_k \left|H_{jk}(\omega) \right|
+% \end{equation}
+
+% The result is shown on figure [[fig:composite_response_function]].
+
+
+figure;
+hold on;
+plot(freqs, squeeze(sum(sum(abs(FRFs)))), '-k');
+plot(freqs, squeeze(sum(sum(abs(FRFs_O)))), '--k');
+hold off;
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([1, 200]);
+
+% Importation of the modal parameters on Matlab
+% Then we import the obtained =.txt= files on Matlab using =readtable= function.
+
+
+shapes_m = readtable('mat/mode_shapes.txt',              'ReadVariableNames', false);  % [Sign / Real / Imag]
+freqs_m  = table2array(readtable('mat/mode_freqs.txt',   'ReadVariableNames', false)); % in [Hz]
+damps_m  = table2array(readtable('mat/mode_damps.txt',   'ReadVariableNames', false)); % in [%]
+modal_a  = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
+modal_b  = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]
+
+
+
+% We guess the number of modes identified from the length of the imported data.
+
+acc_n = 23; % Number of accelerometers
+dir_n = 3; % Number of directions
+dirs = 'XYZ';
+
+mod_n = size(shapes_m,1)/acc_n/dir_n; % Number of modes
+
+
+
+% As the mode shapes are split into 3 parts (direction plus sign, real part and imaginary part), we aggregate them into one array of complex numbers.
+
+T_sign = table2array(shapes_m(:, 1));
+T_real = table2array(shapes_m(:, 2));
+T_imag = table2array(shapes_m(:, 3));
+
+mode_shapes = zeros(mod_n, dir_n, acc_n);
+
+for mod_i = 1:mod_n
+  for acc_i = 1:acc_n
+    % Get the correct section of the signs
+    T = T_sign(acc_n*dir_n*(mod_i-1)+1:acc_n*dir_n*mod_i);
+    for dir_i = 1:dir_n
+      % Get the line corresponding to the sensor
+      i = find(contains(T, sprintf('%i%s',acc_i, dirs(dir_i))), 1, 'first')+acc_n*dir_n*(mod_i-1);
+      mode_shapes(mod_i, dir_i, acc_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
+    end
+  end
+end
+
+% Modal Matrices
+% We would like to arrange the obtained modal parameters into two modal matrices:
+% \[ \Lambda = \begin{bmatrix}
+%   s_1 &        & 0 \\
+%       & \ddots &   \\
+%   0   &        & s_N
+% \end{bmatrix}_{N \times N}; \quad \Psi = \begin{bmatrix}
+%   & & \\
+%   \{\psi_1\} & \dots & \{\psi_N\} \\
+%   & &
+% \end{bmatrix}_{M \times N} \]
+% \[ \{\psi_i\} = \begin{Bmatrix} \psi_{i, 1_x} & \psi_{i, 1_y} & \psi_{i, 1_z} & \psi_{i, 2_x} & \dots & \psi_{i, 23_z} \end{Bmatrix}^T \]
+
+% $M$ is the number of DoF: here it is $23 \times 3 = 69$.
+% $N$ is the number of mode
+
+
+eigen_val_M = diag(2*pi*freqs_m.*(-damps_m/100 + j*sqrt(1 - (damps_m/100).^2)));
+eigen_vec_M = reshape(mode_shapes, [mod_n, acc_n*dir_n]).';
+
+
+
+% Each eigen vector is normalized: $\| \{\psi_i\} \|_2 = 1$
+
+
+% However, the eigen values and eigen vectors appears as complex conjugates:
+% \[ s_r, s_r^*, \{\psi\}_r, \{\psi\}_r^*, \quad r = 1, N \]
+
+% In the end, they are $2N$ eigen values.
+% We then build two extended eigen matrices as follow:
+% \[ \mathcal{S} = \begin{bmatrix}
+%   s_1 &        &     &       &        &      \\
+%       & \ddots &     &       & 0      &      \\
+%       &        & s_N &       &        &      \\
+%       &        &     & s_1^* &        &      \\
+%       & 0      &     &       & \ddots &      \\
+%       &        &     &       &        & s_N^*
+% \end{bmatrix}_{2N \times 2N}; \quad \Phi = \begin{bmatrix}
+%   & & & & &\\
+%   \{\psi_1\} & \dots & \{\psi_N\} & \{\psi_1^*\} & \dots & \{\psi_N^*\} \\
+%   & & & & &
+% \end{bmatrix}_{M \times 2N} \]
+
+
+eigen_val_ext_M = blkdiag(eigen_val_M, conj(eigen_val_M));
+eigen_vec_ext_M = [eigen_vec_M, conj(eigen_vec_M)];
+
+
+
+% We also build the Modal A and Modal B matrices:
+% \begin{equation}
+%   A = \begin{bmatrix}
+%     a_1 &        & 0   \\
+%         & \ddots &     \\
+%     0   &        & a_N
+%   \end{bmatrix}_{N \times N}; \quad B = \begin{bmatrix}
+%     b_1 &        & 0   \\
+%         & \ddots &     \\
+%     0   &        & b_N
+%   \end{bmatrix}_{N \times N}
+% \end{equation}
+% With $a_i$ is the "Modal A" parameter linked to mode i.
+
+
+modal_a_M = diag(complex(modal_a(:, 1), modal_a(:, 2)));
+modal_b_M = diag(complex(modal_b(:, 1), modal_b(:, 2)));
+
+modal_a_ext_M = blkdiag(modal_a_M, conj(modal_a_M));
+modal_b_ext_M = blkdiag(modal_b_M, conj(modal_b_M));
+
+% Matlab Implementation
+
+Hsyn = zeros(69, 69, 801);
+
+for i = 1:length(freqs)
+  Hsyn(:, :, i) = eigen_vec_ext_M*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M.';
+end
+
+
+
+% Because the synthesize frequency response functions are representing the displacement response in $[m/N]$, we multiply each element of the FRF matrix by $(j \omega)^2$ in order to obtain the acceleration response in $[m/s^2/N]$.
+
+
+for i = 1:size(Hsyn, 1)
+  Hsyn(i, :, :) = squeeze(-Hsyn(i, :, :)).*(j*2*pi*freqs).^2;
+end
+
+% Original and Synthesize FRF matrix comparison
+
+acc_o = 15; dir_o = 1; dir_i = 1;
+
+figure;
+ax1 = subplot(2, 1, 1);
+hold on;
+plot(freqs, abs(squeeze(FRFs(3*(acc_o-1)+dir_o,          dir_i, :))), 'DisplayName', 'Original');
+plot(freqs, abs(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :))), 'DisplayName', 'Synthesize');
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+legend('location', 'northwest');
+
+ax2 = subplot(2, 1, 2);
+hold on;
+plot(freqs, mod(180/pi*phase(squeeze(FRFs(3*(acc_o-1)+dir_o,          dir_i, :)))+180, 360)-180);
+plot(freqs, mod(180/pi*phase(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :)))+180, 360)-180);
+hold off;
+yticks(-360:90:360); ylim([-180, 180]);
+set(gca, 'xscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+linkaxes([ax1,ax2],'x');
+xlim([1, 200]);
+
+% Synthesize FRF that has not yet been measured
+
+dir_names = {'X', 'Y', 'Z'};
+
+accs = [1]; dirs = [1:3];
+
+figure;
+ax1 = subplot(2, 1, 1);
+hold on;
+for acc_i = accs
+  for dir_i = dirs
+    plot(freqs, abs((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :))), 'DisplayName', sprintf('Acc %i - %s', acc_i, dir_names{dir_i}));
+  end
+end
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('Magnitude [$\frac{m}{N}$]');
+legend('location', 'southwest');
+
+ax2 = subplot(2, 1, 2);
+hold on;
+for acc_i = accs
+  for dir_i = dirs
+    plot(freqs, mod(180/pi*phase((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :)))+180, 360)-180);
+  end
+end
+hold off;
+yticks(-360:90:360); ylim([-180, 180]);
+set(gca, 'xscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+linkaxes([ax1,ax2],'x');
+xlim([1, 200]);
+
+% Modal Complexity
+%   <<sec:modal_complexity>>
+% A method of displaying *modal complexity* is by plotting the elements of the eigenvector on an *Argand diagram* (complex plane), such as the ones shown in figure [[fig:modal_complexity_small]].
+
+% To evaluate the complexity of the modes, we plot a polygon around the extremities of the individual vectors.
+% The obtained area of this polygon is then compared with the area of the circle which is based on the length of the largest vector element. The resulting ratio is used as an *indication of the complexity of the mode*.
+
+% A mode with small complexity is shown on figure [[fig:modal_complexity_small]] whereas an highly complex mode is shown on figure [[fig:modal_complexity_high]].
+% The complexity of all the modes are compared on figure [[fig:modal_complexities]].
+
+
+figure;
+mod_i = 1;
+i_max = convhull(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)));
+radius = max(abs(eigen_vector_M(:, mod_i)));
+theta = linspace(0, 2*pi, 100);
+
+hold on;
+plot(radius*cos(theta), radius*sin(theta), '-');
+plot(real(eigen_vector_M(i_max, mod_i)), imag(eigen_vector_M(i_max, mod_i)), '-');
+plot(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)), 'ko');
+hold off;
+xlabel('Real Part'); ylabel('Imaginary Part');
+title(sprintf('Mode %i', mod_i));
+axis manual equal
+
+
+
+% #+NAME: fig:modal_complexity_small
+% #+CAPTION: Modal Complexity of one mode with small complexity
+% [[file:figs/modal_complexity_small.png]]
+
+
+mod_i = 8;
+i_max = convhull(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)));
+radius = max(abs(eigen_vector_M(:, mod_i)));
+theta = linspace(0, 2*pi, 100);
+
+figure;
+hold on;
+plot(radius*cos(theta), radius*sin(theta), '-');
+plot(real(eigen_vector_M(i_max, mod_i)), imag(eigen_vector_M(i_max, mod_i)), '-');
+plot(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)), 'ko');
+hold off;
+xlabel('Real Part'); ylabel('Imaginary Part');
+title(sprintf('Mode %i', mod_i));
+axis manual equal
+
+
+
+% #+NAME: fig:modal_complexity_high
+% #+CAPTION: Modal Complexity of one higly complex mode
+% [[file:figs/modal_complexity_high.png]]
+
+
+modes_complexity = zeros(mod_n, 1);
+for mod_i = 1:mod_n
+  i = convhull(real(eigen_vector_M(:, mod_i)), imag(eigen_vector_M(:, mod_i)));
+  area_complex = polyarea(real(eigen_vector_M(i, mod_i)), imag(eigen_vector_M(i, mod_i)));
+  area_circle = pi*max(abs(eigen_vector_M(:, mod_i)))^2;
+  modes_complexity(mod_i) = area_complex/area_circle;
+end
+
+figure;
+plot(1:mod_n, modes_complexity, 'ok');
+ylim([0, 1]);
+xlabel('Mode Number'); ylabel('Modal Complexity');
+
+% Modal Assurance Criterion (MAC)
+% The MAC is calculated as the normalized scalar product of the two sets of vectors $\{\psi_A\}$ and $\{\psi_X\}$.
+% The resulting scalars are arranged into the MAC matrix cite:pastor12_modal_assur_criter:
+% \begin{equation}
+%   \text{MAC}(r, q) = \frac{\left| \{\psi_A\}_r^T \{\psi_X\}_q^* \right|^2}{\left( \{\psi_A\}_r^T \{\psi_A\}_r^* \right) \left( \{\psi_X\}_q^T \{\psi_X\}_q^* \right)}
+% \end{equation}
+
+% An equivalent formulation is:
+% \begin{equation}
+%   \text{MAC}(r, q) = \frac{\left| \sum_{j=1}^n \{\psi_A\}_j \{\psi_X\}_j^* \right|^2}{\left( \sum_{j=1}^n |\{\psi_A\}_j|^2 \right) \left( \sum_{j=1}^n |\{\psi_X\}_j|^2 \right)}
+% \end{equation}
+
+% The MAC takes value between 0 (representing no consistent correspondence) and 1 (representing a consistent correspondence).
+
+% We compute the autoMAC matrix that compares all the possible combinations of mode shape pairs for only one set of mode shapes. The result is shown on figure [[fig:automac]].
+
+
+autoMAC = eye(size(eigen_vec_M, 2));
+
+for r = 1:size(eigen_vec_M, 2)
+  for q = r+1:size(eigen_vec_M, 2)
+    autoMAC(r, q) = abs(eigen_vec_M(r, :)*eigen_vec_M(q, :)')^2/((eigen_vec_M(r, :)*eigen_vec_M(r, :)')*(eigen_vec_M(q, :)*eigen_vec_M(q, :)'));
+    autoMAC(q, r) = autoMAC(r, q);
+  end
+end
+
+figure;
+imagesc(autoMAC);
+colormap('parula');
+colorbar;
+xticks(1:1:size(eigen_val_M, 1)); yticks(1:1:size(eigen_val_M, 1));
+xticklabels(num2str(round(imag(diag(eigen_val_M))/2/pi, 1)));
+xtickangle(90);
+yticklabels(num2str(round(imag(diag(eigen_val_M))/2/pi, 1)));
+set(gca,'YDir','normal');
+
+% Matlab Implementation
+% The obtained mode shapes matrix that gives the mode shapes of each solid bodies with respect to the fixed frame $\{O\}$, =mode_shapes_O=, is an $n \times p \times q$ with:
+% - $n$ is the number of modes: 21
+% - $p$ is the number of DOFs for each solid body: 6
+% - $q$ is the number of solid bodies: 6
+
+
+mode_shapes_O = zeros(mod_n, 6, length(solid_names));
+
+for mod_i = 1:mod_n
+  for solid_i = 1:length(solid_names)
+    acc_i = solids.(solid_names{solid_i});
+
+    Y = mode_shapes(mod_i, :, acc_i);
+    Y = Y(:);
+
+    A = zeros(3*length(acc_i), 6);
+    for i = 1:length(acc_i)
+      A(3*(i-1)+1:3*i, 1:3) = eye(3);
+
+      A(3*(i-1)+1:3*i, 4:6) = [ 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];
+    end
+
+    mode_shapes_O(mod_i, :, solid_i) = A\Y;
+  end
+end
+
+
+
+% We then rearrange the eigen vectors in another way:
+% \[ \Psi_O = \begin{bmatrix}
+%   & & \\
+%   \{\psi_1\} & \dots & \{\psi_N\} \\
+%   & &
+% \end{bmatrix}_{M \times N} \]
+% with
+% \[ \{\psi\}_r = \begin{Bmatrix}
+%   \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z}
+% \end{Bmatrix}^T \]
+
+% With $M = 6 \times 6$ is the new number of DOFs and $N=21$ is the number of modes.
+
+
+eigen_vec_O = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)]).';
+
+eigen_vec_ext_O = [eigen_vec_O, conj(eigen_vec_O)];
+
+% Matlab Implementation
+
+Hsyn_O = zeros(36, 36, 801);
+
+for i = 1:length(freqs)
+  Hsyn_O(:, :, i) = eigen_vec_ext_O*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_O.';
+end
+
+
+
+% Because the synthesize frequency response functions are representing the displacement response in $[m/N]$, we multiply each element of the FRF matrix by $(j \omega)^2$ in order to obtain the acceleration response in $[m/s^2/N]$.
+
+
+for i = 1:size(Hsyn_O, 1)
+  Hsyn_O(i, :, :) = -squeeze(Hsyn_O(i, :, :)).*(j*2*pi*freqs).^2;
+end
+
+% TODO Test to have the original inputs
+
+Hsyn_O = zeros(36, 3, 801);
+
+for i = 1:length(freqs)
+  Hsyn_O(:, :, i) = eigen_vec_ext_O*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M(10*3+1:11*3, :).';
+end
+
+for i = 1:size(Hsyn_O, 1)
+  Hsyn_O(i, :, :) = -squeeze(Hsyn_O(i, :, :)).*(j*2*pi*freqs).^2;
+end
+
+% TODO Verification
+
+solid_o = 6; dir_o = 3;
+solid_i = 3; dir_i = 3;
+
+figure;
+ax1 = subplot(2, 1, 1);
+hold on;
+plot(freqs, abs(squeeze(FRFs_O((solid_o-1)*6+dir_o, dir_i, :))), 'DisplayName', 'Original');
+plot(freqs, abs(squeeze(Hsyn_O((solid_o-1)*6+dir_o, dir_i, :))), 'DisplayName', 'Synthesize');
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+legend('location', 'northwest');
+
+ax2 = subplot(2, 1, 2);
+hold on;
+plot(freqs, mod(180/pi*phase(squeeze(FRFs_O((solid_o-1)*6+dir_o, dir_i, :)))+180, 360)-180);
+plot(freqs, mod(180/pi*phase(squeeze(Hsyn_O((solid_o-1)*6+dir_o, dir_i, :)))+180, 360)-180);
+hold off;
+yticks(-360:90:360);
+set(gca, 'xscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+linkaxes([ax1,ax2],'x');
+xlim([1, 200]);
+
+% New FRFs
+
+dof_names = {'X', 'Y', 'Z', '\theta_X', '\theta_Y', '\theta_Z'};
+solid_i = 6; dirs = [1:3];
+
+figure;
+ax1 = subplot(2, 1, 1);
+hold on;
+for dir_i = dirs
+  plot(freqs, abs(squeeze(Hsyn_O((solid_i-1)*6+dir_i, (solid_i-1)*6+dir_i, :))), 'DisplayName', sprintf('Solid %s - %s', solid_names{solid_i}, dof_names{dir_i}));
+end
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+legend('location', 'northwest');
+
+ax2 = subplot(2, 1, 2);
+hold on;
+for dir_i = dirs
+  plot(freqs, mod(180/pi*phase(squeeze(Hsyn_O((solid_i-1)*6+dir_i, (solid_i-1)*6+dir_i, :)))+180, 360)-180);
+end
+hold off;
+yticks(-360:90:360);
+set(gca, 'xscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+linkaxes([ax1,ax2],'x');
+xlim([1, 200]);
+
+% Compare Mode Shapes
+% Let's say we want to see for the first mode which DOFs can be neglected.
+% In order to do so, we should estimate the motion of each stage in particular directions.
+% If we look at the z motion for instance, we will find that we cannot neglect that motion (because of the tilt causing z motion).
+
+
+mode_i = 3;
+dof_i = 6;
+
+mode = eigen_vector_M(dof_i:6:end, mode_i);
+
+figure;
+hold on;
+for i=1:length(mode)
+  plot([0, real(mode(i))], [0, imag(mode(i))], '-', 'DisplayName', solid_names{i});
+end
+hold off;
+legend();
+
+figure;
+subplot(2, 1, 1);
+hold on;
+for i=1:length(mode)
+  plot(1, norm(mode(i)), 'o');
+end
+hold off;
+ylabel('Amplitude');
+
+subplot(2, 1, 2);
+hold on;
+for i=1:length(mode)
+  plot(1, 180/pi*angle(mode(i)), 'o', 'DisplayName', solid_names{i});
+end
+hold off;
+ylim([-180, 180]); yticks([-180:90:180]);
+ylabel('Phase [deg]');
+legend();
+
+test = mode_shapes_O(10, 1, :)/norm(squeeze(mode_shapes_O(10, 1, :)));
+test = mode_shapes_O(10, 2, :)/norm(squeeze(mode_shapes_O(10, 2, :)));
+
+% Modes shapes transformation
+% First, we load the position of the Center of Mass of each solid body with respect to the point of interest.
+
+model_com = table2array(readtable('mat/model_solidworks_com.txt', 'ReadVariableNames', false));
+
+
+
+% Then, we reshape the data.
+
+model_com = reshape(model_com , [3, 6]);
+
+
+
+% Now we convert the mode shapes expressed in the DOF of the accelerometers to the DoF of each solid body (translations and rotations) with respect to their own CoM.
+
+mode_shapes_CoM = zeros(mod_n, 6, length(solid_names));
+
+for mod_i = 1:mod_n
+  for solid_i = 1:length(solid_names)
+    acc_i = solids.(solid_names{solid_i});
+
+    Y = mode_shapes(mod_i, :, acc_i);
+    Y = Y(:);
+
+    A = zeros(3*length(acc_i), 6);
+    for i = 1:length(acc_i)
+      acc_pos_com = acc_pos(i, :).' - model_com(:, i);
+
+      A(3*(i-1)+1:3*i, 1:3) = eye(3);
+      A(3*(i-1)+1:3*i, 4:6) = [ 0               acc_pos_com(3) -acc_pos_com(2) ;
+                               -acc_pos_com(3)  0               acc_pos_com(1) ;
+                                acc_pos_com(2) -acc_pos_com(1)  0];
+    end
+
+    mode_shapes_CoM(mod_i, :, solid_i) = A\Y;
+  end
+end
+
+
+
+% We then rearrange the eigen vectors in another way:
+% \[ \Psi_{\text{CoM}} = \begin{bmatrix}
+%   & & \\
+%   \{\psi_1\} & \dots & \{\psi_N\} \\
+%   & &
+% \end{bmatrix}_{M \times N} \]
+% with
+% \[ \{\psi\}_r = \begin{Bmatrix}
+%   \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z}
+% \end{Bmatrix}^T \]
+
+% With $M = 6 \times 6$ is the new number of DOFs and $N=21$ is the number of modes.
+
+% Each eigen vector is normalized.
+
+
+eigen_vec_CoM = reshape(mode_shapes_CoM, [mod_n, 6*length(solid_names)]).';
+% eigen_vec_CoM = eigen_vec_CoM./vecnorm(eigen_vec_CoM);
+eigen_vec_ext_CoM = [eigen_vec_CoM, conj(eigen_vec_CoM)];
+
+model_mass = table2array(readtable('mat/model_solidworks_mass.txt', 'ReadVariableNames', false));
+model_inertia = table2array(readtable('mat/model_solidworks_inertia.txt', 'ReadVariableNames', false));
+
+
+
+% \[ M = \begin{bmatrix}
+%   M_1 &     &     &     &     &     \\
+%       & M_2 &     &     &     &     \\
+%       &     & M_3 &     &     &     \\
+%       &     &     & M_4 &     &     \\
+%       &     &     &     & M_5 &     \\
+%       &     &     &     &     & M_6 \\
+% \end{bmatrix}, \text{ with } M_i = \begin{bmatrix}
+%   m_i &     &     &          &          &          \\
+%       & m_i &     &          &          &          \\
+%       &     & m_i &          &          &          \\
+%       &     &     & I_{i,xx} & I_{i,yx} & I_{i,zx} \\
+%       &     &     & I_{i,xy} & I_{i,yy} & I_{i,zy} \\
+%       &     &     & I_{i,xz} & I_{i,yz} & I_{i,zz} \\
+% \end{bmatrix} \]
+
+
+M = zeros(6*6, 6*6);
+
+for i = 1:6
+  M((i-1)*6+1:i*6, (i-1)*6+1:i*6) = blkdiag(model_mass(i)*eye(3), model_inertia((i-1)*3+1:i*3, 1:3));
+end
+
+% Mass-normalized Eigen Vectors
+% To do so, it seems that we need to know the mass matrix $[M]$. Then:
+% \[ \{\phi\}_r = \frac{1}{\sqrt{m_r}} \{\psi\}_r \text{ where } m_r = \{\psi\}_r^T [M] \{\psi\}_r \]
+
+% Mass-normalized eigen vectors are very important for the synthesis and spatial model extraction.
+
+% Mass Matrix could be estimated from the SolidWorks model.
+
+% We compute the modal masses that will be used for normalization.
+
+mr = zeros(size(eigen_vector_matrix_CoM, 2), 1);
+
+for i = 1:length(mr)
+  mr(i) = real(eigen_vector_matrix_CoM(:, i).' * M * eigen_vector_matrix_CoM(:, i));
+end
+
+figure;
+plot(freqs_m, mr, 'ko');
+xlabel('Mode Frequency [Hz]'); ylabel('Modal mass [kg]');
+
+
+
+% #+NAME: fig:modal_mass
+% #+CAPTION: Modal masses
+% [[file:figs/modal_mass.png]]
+
+% And finally, we compute the mass-normalized eigen vectors.
+
+eigen_vector_mass_CoM = zeros(size(eigen_vector_matrix_CoM));
+
+for i = 1:size(eigen_vector_matrix_CoM, 2)
+  eigen_vector_mass_CoM(:, i) = 1/sqrt(mr(i)) * eigen_vector_matrix_CoM(:, i);
+end
+
+
+
+% Verification that
+% \[ [\Phi]^T [M] [\Phi] = [I] \]
+
+
+eigen_vector_matrix_CoM.'*M*eigen_vector_matrix_CoM
+
+eigen_vector_mass_CoM*M*eigen_vector_mass_CoM.'
+
+
+
+% Other test for normalized eigen vectors
+
+eigen_vector_mass_CoM = (eigen_vector_matrix_CoM.'*diag(diag(M))*eigen_vector_matrix_CoM)^(-0.5) * eigen_vector_matrix_CoM';
+eigen_vector_mass_CoM = eigen_vector_mass_CoM.';
+
+% Full Response Model from modal model (synthesis)
+% In general, the form of response model with which we are concerned is an *FRF matrix* whose order is dictated by the number of coordinates $n$ included in the test.
+% Also, as explained, it is normal in practice to measured and to analyze just a *subset of the FRF matrix* but rather to measure the full FRF matrix.
+% Usually *one column* or *one row* with a few additional elements are measured.
+
+% Thus, if we are to construct an acceptable response model, it will be necessary to synthesize those elements which have not been directly measured.
+% However, in principle, this need present no major problem as it is possible to compute the full FRF matrix from a modal model using:
+% \begin{equation}
+%   [H]_{n\times n} = [\Phi]_{n\times m} [\lambda_r^2 - \omega^2]_{m\times m}^{-1} [\Phi]_{m\times n}^T
+% \end{equation}
+
+% $\{\Phi\}$ is a *mass-normalized* eigen vector.
+
+
+FRF_matrix_CoM = zeros(size(eigen_vector_mass_CoM, 1), size(eigen_vector_mass_CoM, 1), length(freqs));
+
+for i = 1:length(freqs)
+  FRF_matrix_CoM(:, :, i) = eigen_vector_mass_CoM*inv(eigen_value_M^2 - (2*pi*freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass_CoM.';
+end
+
+exc_dir = 3;
+meas_mass = 6;
+meas_dir = 3;
+
+figure;
+hold on;
+plot(freqs, abs(squeeze(FRF_matrix_CoM((meas_mass-1)*6 + meas_dir, 6*2+exc_dir, :))));
+plot(freqs, abs(squeeze(FRFs_CoM((meas_mass-1)*6 + meas_dir, exc_dir, :))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([1, 200]);
+
+FRF_matrix = zeros(size(eigen_vector_mass, 1), size(eigen_vector_mass, 1), length(freqs));
+
+for i = 1:length(freqs)
+  FRF_matrix(:, :, i) = eigen_vector_mass*inv(eigen_value_M - (freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass.';
+end
+
+figure;
+hold on;
+plot(freqs, abs(squeeze(FRFs_O(1, 1, :))));
+plot(freqs, abs(squeeze(FRF_matrix(1, 1, :))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+xlim([1, 200]);

modal-analysis/matlab/modal_frf_coh.m

@@ -224,7 +224,7 @@ ax1 = subplot(2, 1, 1);
 plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_AvXSpc_2_1_RMS_Y_Mod);
 set(gca, 'Yscale', 'log');
 set(gca, 'XTickLabel',[]);
-ylabel('Magnitude');
+ylabel('Magnitude [$\frac{m/s^2}{N}$]');
 
 ax2 = subplot(2, 1, 2);
 plot(meas1.FFT1_AvSpc_2_RMS_X_Val, meas1.FFT1_AvXSpc_2_1_RMS_Y_Phas);
@@ -255,7 +255,7 @@ ylabel('Coherence');
 % - $q$ is the number of frequency points $\omega_i$
 
 % Thus, the FRF matrix is an $69 \times 3 \times 801$ matrix.
-
+% We do the same thing for the coherence matrix.
 
 % #+begin_important
 % For each frequency point $\omega_i$, we obtain a 2D matrix:
@@ -372,23 +372,23 @@ save('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
 % #+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 motion of the rigid body of figure [[fig:aligned_accelerometers]] is defined by its displacement $\delta p$ 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
+%   \delta p_1 &= \delta p + \Omega \times p_1 \\
+%   \delta p_2 &= \delta p + \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}
+%   \delta p_{x1} &= \delta p_x + \Omega_y p_{z1} - \Omega_z p_{y1} \\
+%   \delta p_{x2} &= \delta p_x + \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}
+%   \delta p_{x1} = \delta p_{x2}
 % \end{equation}
 
 % #+begin_important
@@ -430,7 +430,7 @@ for i = 1:size(acc_i, 1)
   end
 
   if rem(i, 3) == 1
-    ylabel('Amplitude');
+    ylabel('Amplitude [$\frac{m/s^2}{N}$]');
   end
   xlim([1, 200]);
   title(sprintf('Acc %i and %i - X', acc_i(i, 1), acc_i(i, 2)));
@@ -472,8 +472,44 @@ for i = 1:size(acc_i, 1)
   end
 
   if rem(i, 3) == 1
-    ylabel('Amplitude');
+    ylabel('Amplitude [$\frac{m/s^2}{N}$]');
   end
   xlim([1, 200]);
   title(sprintf('Acc %i and %i - Y', acc_i(i, 1), acc_i(i, 2)));
 end
+
+% Verification of the principle of reciprocity
+% Because we expect our system to follow the principle of reciprocity.
+% That is to say the response $X_j$ at some degree of freedom $j$ due to a force $F_k$ applied on DOF $k$ should be the same as the response $X_k$ due to a force $F_j$:
+% \[ H_{jk} = \frac{X_j}{F_k} = \frac{X_k}{F_j} = H_{kj} \]
+
+% This comes from the fact that we expect to have symmetric mass, stiffness and damping matrices.
+
+% In order to access the quality of the data and the validity of the measured FRF, we then check that the reciprocity between $H_{jk}$ and $H_{kj}$ is of an acceptable level.
+
+% We can verify this reciprocity using 3 different pairs of response/force.
+
+
+dir_names = {'X', 'Y', 'Z'};
+
+figure;
+for i = 1:3
+  subplot(3, 1, i)
+  a = mod(i, 3)+1;
+  b = mod(i-2, 3)+1;
+  hold on;
+  plot(freqs, abs(squeeze(FRFs(3*(11-1)+a, b, :))), 'DisplayName', sprintf('$\\frac{F_%s}{D_%s}$', dir_names{a}, dir_names{b}));
+  plot(freqs, abs(squeeze(FRFs(3*(11-1)+b, a, :))), 'DisplayName', sprintf('$\\frac{F_%s}{D_%s}$', dir_names{b}, dir_names{a}));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  if i == 3
+    xlabel('Frequency [Hz]');
+  else
+    set(gca, 'XTickLabel',[]);
+  end
+  if i == 2
+    ylabel('Amplitude [$\frac{m/s^2}{N}$]');
+  end
+  xlim([1, 200]);
+  legend('location', 'northwest');
+end

modal-analysis/modal_extraction.org

@@ -43,8 +43,7 @@
 
 * Introduction                                                       :ignore:
 The goal here is to extract the modal parameters describing the modes of station being studied, namely:
-- the eigen frequencies
-- the modal damping
+- the eigen frequencies and the modal damping (eigen values)
 - the mode shapes (eigen vectors)
 This is done from the FRF matrix previously extracted from the measurements.
 
@@ -61,10 +60,12 @@ However, there are multiple level of complexity:
 - works on a set of many FRF plots all obtained from the same structure
 The third method is the most complex but gives better results. This is the one we will use in section [[sec:modal_extraction]].
 
-
 From the modal model, it is possible to obtain a graphic display of the mode shapes (section [[sec:mode_shape_display]]).
 
 
+In order to validate the quality of the modal model, we will synthesize the FRF matrix from the modal model and compare it with the FRF measured (section [[sec:validity_measurement]]).
+
+
 The modes of the structure are expected to be complex, however real modes are easier to work with when it comes to obtain a spatial model from the modal parameters.
 We will thus study the complexity of those modes, in section [[sec:modal_complexity]], and see if we can estimate real modes from the complex modes.
 
@@ -73,6 +74,8 @@ The mode obtained from the modal software describe the motion of the structure a
 However, we would like to describe the motion of each stage (solid body) of the structure in its 6 DOFs.
 This in done in section [[sec:global_mode_shapes]].
 
+Residues will have to be included.
+
 * ZIP file containing the data and matlab files                      :ignore:
 #+begin_src bash :exports none :results none
   if [ matlab/modal_extraction.m -nt data/modal_extraction.zip ]; then
@@ -106,6 +109,7 @@ This in done in section [[sec:global_mode_shapes]].
 * Load Data
 #+begin_src matlab
   load('mat/frf_coh_matrices.mat', 'FRFs', 'COHs', 'freqs');
+  load('mat/geometry.mat', 'solids', 'solid_names', 'acc_pos');
   load('mat/frf_o.mat', 'FRFs_O');
 #+end_src
 
@@ -240,13 +244,13 @@ The result is shown on figure [[fig:composite_response_function]].
 * Modal parameter extraction
   <<sec:modal_extraction>>
 ** OROS - Modal software
-Modal identification can be done within the [[https://www.oros.com/solutions/structural-dynamics/modal-analysis/][Modal software of OROS]]. The manual for the software is available [[file:data/oros_modal_manual.pdf][here]].
+Modal identification are done within the [[https://www.oros.com/solutions/structural-dynamics/modal-analysis/][Modal software of OROS]]. The manual for the software is available [[file:data/oros_modal_manual.pdf][here]].
 
 Several modal parameter extraction methods are available.
-We choose to use the "broad band method" as it permits to identify the modal parameters using all the FRF curves at the same times.
+We choose to use the "*broad band*" method as it permits to identify the modal parameters using all the FRF curves at the same time.
 It takes into account the fact the the properties of all the individual curves are related by being from the same structure: all FRF plots on a given structure should indicate the same values for the natural frequencies and damping factor of each mode.
 
-Such method also have the advantage of producing a unique and consistent model as direct output.
+Such method also have the advantage of producing a *unique and consistent model* as direct output.
 
 In order to apply this method, we select the frequency range of interest and we give an estimate of how many modes are present.
 
@@ -254,16 +258,16 @@ Then, it shows a stabilization charts, such as the one shown on figure [[fig:sta
 
 #+name: fig:stabilization_chart
 #+caption: Stabilization Chart
-[[file:figs/stabilisation_chart.jpg]]
+[[file:img/modal_software/stabilisation_chart.jpg]]
 
 We can then run the modal analysis, and the software will identify the modal damping and mode shapes at the selected frequency modes.
 
 ** Exported modal parameters
 The obtained modal parameters are:
-- the resonance frequencies
-- the modes shapes
-- the modal damping
-- the residues
+- resonance frequencies in Hertz
+- modal damping ratio in percentage
+- (complex) modes shapes for each measured DoF
+- modal A and modal B which are parameters important for further normalization
 
 They are all exported in a text file named =modes.asc=.
 Its first 20 lines as shown below.
@@ -308,25 +312,23 @@ We split this big =modes.asc= file into sub text files using =bash=. The obtaine
 
 #+name: tab:modes_files
 #+caption: Split =modes.asc= file
-| Filename               | Content                                          |
-|------------------------+--------------------------------------------------|
-| =mat/mode_shapes.txt=  | mode shapes                                      |
-| =mat/mode_freqs.txt=   | resonance frequencies                            |
-| =mat/mode_damps.txt=   | modal damping                                    |
-| =mat/mode_modal_a.txt= | modal residues at low frequency (to be checked)  |
-| =mat/mode_modal_b.txt= | modal residues at high frequency (to be checked) |
+| Filename               | Content                    |
+|------------------------+----------------------------|
+| =mat/mode_shapes.txt=  | mode shapes                |
+| =mat/mode_freqs.txt=   | resonance frequencies [Hz] |
+| =mat/mode_damps.txt=   | modal damping [per]        |
+| =mat/mode_modal_a.txt= | Modal A                    |
+| =mat/mode_modal_b.txt= | Modal B                    |
 
 ** Importation of the modal parameters on Matlab
 Then we import the obtained =.txt= files on Matlab using =readtable= function.
-#+begin_src matlab
-  shapes  = readtable('mat/mode_shapes.txt',              'ReadVariableNames', false);  % [Sign / Real / Imag]
-  freqs   = table2array(readtable('mat/mode_freqs.txt',   'ReadVariableNames', false)); % in [Hz]
-  damps   = table2array(readtable('mat/mode_damps.txt',   'ReadVariableNames', false)); % in [%]
-  modal_a = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
-  modal_b = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]
 
-  modal_a = complex(modal_a(:, 1), modal_a(:, 2));
-  modal_b = complex(modal_b(:, 1), modal_b(:, 2));
+#+begin_src matlab
+  shapes_m = readtable('mat/mode_shapes.txt',              'ReadVariableNames', false);  % [Sign / Real / Imag]
+  freqs_m  = table2array(readtable('mat/mode_freqs.txt',   'ReadVariableNames', false)); % in [Hz]
+  damps_m  = table2array(readtable('mat/mode_damps.txt',   'ReadVariableNames', false)); % in [%]
+  modal_a  = table2array(readtable('mat/mode_modal_a.txt', 'ReadVariableNames', false)); % [Real / Imag]
+  modal_b  = table2array(readtable('mat/mode_modal_b.txt', 'ReadVariableNames', false)); % [Real / Imag]
 #+end_src
 
 We guess the number of modes identified from the length of the imported data.
@@ -335,16 +337,16 @@ We guess the number of modes identified from the length of the imported data.
   dir_n = 3; % Number of directions
   dirs = 'XYZ';
 
-  mod_n = size(shapes,1)/acc_n/dir_n; % Number of modes
+  mod_n = size(shapes_m,1)/acc_n/dir_n; % Number of modes
 #+end_src
 
 As the mode shapes are split into 3 parts (direction plus sign, real part and imaginary part), we aggregate them into one array of complex numbers.
 #+begin_src matlab
-  T_sign = table2array(shapes(:, 1));
-  T_real = table2array(shapes(:, 2));
-  T_imag = table2array(shapes(:, 3));
+  T_sign = table2array(shapes_m(:, 1));
+  T_real = table2array(shapes_m(:, 2));
+  T_imag = table2array(shapes_m(:, 3));
 
-  modes = zeros(mod_n, acc_n, dir_n);
+  mode_shapes = zeros(mod_n, dir_n, acc_n);
 
   for mod_i = 1:mod_n
     for acc_i = 1:acc_n
@@ -353,7 +355,7 @@ As the mode shapes are split into 3 parts (direction plus sign, real part and im
       for dir_i = 1:dir_n
         % Get the line corresponding to the sensor
         i = find(contains(T, sprintf('%i%s',acc_i, dirs(dir_i))), 1, 'first')+acc_n*dir_n*(mod_i-1);
-        modes(mod_i, acc_i, dir_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
+        mode_shapes(mod_i, dir_i, acc_i) = str2num([T_sign{i}(end-1), '1'])*complex(T_real(i),T_imag(i));
       end
     end
   end
@@ -361,7 +363,7 @@ As the mode shapes are split into 3 parts (direction plus sign, real part and im
 
 The obtained mode frequencies and damping are shown below.
 #+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([(1:length(freqs))', freqs, damps], {}, {'Mode number', 'Frequency [Hz]', 'Damping [%]'}, ' %.1f ');
+  data2orgtable([(1:length(freqs_m))', freqs_m, damps_m], {}, {'Mode number', 'Frequency [Hz]', 'Damping [%]'}, ' %.1f ');
 #+end_src
 
 #+name: tab:obtained_modes_freqs_damps
@@ -391,26 +393,84 @@ The obtained mode frequencies and damping are shown below.
 |        20.0 |          150.1 |         2.2 |
 |        21.0 |          164.7 |         1.4 |
 
+** Theory
+It seems that the modal analysis software makes the *assumption* of viscous damping for the model with which it tries to fit the FRF measurements.
+
+If we note $N$ the number of modes identified, then there are $2N$ eigenvalues and eigenvectors given by the software:
+\begin{align}
+  s_r &= \omega_r (-\xi_r + i \sqrt{1 - \xi_r^2}),\quad s_r^* \\
+  \{\psi_r\} &= \begin{Bmatrix} \psi_{1_x} & \psi_{2_x} & \dots & \psi_{23_x} & \psi_{1_y} & \dots & \psi_{1_z} & \dots & \psi_{23_z} \end{Bmatrix}^T, \quad \{\psi_r\}^*
+\end{align}
+for $r = 1, \dots, N$ where $\omega_r$ is the natural frequency and $\xi_r$ is the critical damping ratio for that mode.
+
 ** Modal Matrices
-We arrange the obtained modal parameters into matrices:
-\[ \Omega = \begin{bmatrix}
-  \omega_1^2 &        & 0 \\
-             & \ddots &   \\
-  0          &        & \omega_n^2
-\end{bmatrix}; \quad \Psi = \begin{bmatrix}
+We would like to arrange the obtained modal parameters into two modal matrices:
+\[ \Lambda = \begin{bmatrix}
+  s_1 &        & 0 \\
+      & \ddots &   \\
+  0   &        & s_N
+\end{bmatrix}_{N \times N}; \quad \Psi = \begin{bmatrix}
   & & \\
-  \{\psi_1\} & \dots & \{\psi_n\} \\
+  \{\psi_1\} & \dots & \{\psi_N\} \\
   & &
-\end{bmatrix} \]
-with $n$ the number of identified modes and:
-\[ \{\psi_1\} = \begin{Bmatrix} \psi_{1_x} & \psi_{2_x} & \dots & \psi_{23_x} & \psi_{1_y} & \dots & \psi_{1_z} & \dots & \psi_{23_z} \end{Bmatrix}^T \]
+\end{bmatrix}_{M \times N} \]
+\[ \{\psi_i\} = \begin{Bmatrix} \psi_{i, 1_x} & \psi_{i, 1_y} & \psi_{i, 1_z} & \psi_{i, 2_x} & \dots & \psi_{i, 23_z} \end{Bmatrix}^T \]
+
+$M$ is the number of DoF: here it is $23 \times 3 = 69$.
+$N$ is the number of mode
 
 #+begin_src matlab
-  eigen_value_M = diag(freqs*2*pi);
-  eigen_vector_M = reshape(modes, [mod_n, acc_n*dir_n])';
+  eigen_val_M = diag(2*pi*freqs_m.*(-damps_m/100 + j*sqrt(1 - (damps_m/100).^2)));
+  eigen_vec_M = reshape(mode_shapes, [mod_n, acc_n*dir_n]).';
 #+end_src
 
-Each eigen vector is normalized: $\| \{\psi_i\}\} \|_2 = 1$
+Each eigen vector is normalized: $\| \{\psi_i\} \|_2 = 1$
+
+
+However, the eigen values and eigen vectors appears as complex conjugates:
+\[ s_r, s_r^*, \{\psi\}_r, \{\psi\}_r^*, \quad r = 1, N \]
+
+In the end, they are $2N$ eigen values.
+We then build two extended eigen matrices as follow:
+\[ \mathcal{S} = \begin{bmatrix}
+  s_1 &        &     &       &        &      \\
+      & \ddots &     &       & 0      &      \\
+      &        & s_N &       &        &      \\
+      &        &     & s_1^* &        &      \\
+      & 0      &     &       & \ddots &      \\
+      &        &     &       &        & s_N^*
+\end{bmatrix}_{2N \times 2N}; \quad \Phi = \begin{bmatrix}
+  & & & & &\\
+  \{\psi_1\} & \dots & \{\psi_N\} & \{\psi_1^*\} & \dots & \{\psi_N^*\} \\
+  & & & & &
+\end{bmatrix}_{M \times 2N} \]
+
+#+begin_src matlab
+  eigen_val_ext_M = blkdiag(eigen_val_M, conj(eigen_val_M));
+  eigen_vec_ext_M = [eigen_vec_M, conj(eigen_vec_M)];
+#+end_src
+
+We also build the Modal A and Modal B matrices:
+\begin{equation}
+  A = \begin{bmatrix}
+    a_1 &        & 0   \\
+        & \ddots &     \\
+    0   &        & a_N
+  \end{bmatrix}_{N \times N}; \quad B = \begin{bmatrix}
+    b_1 &        & 0   \\
+        & \ddots &     \\
+    0   &        & b_N
+  \end{bmatrix}_{N \times N}
+\end{equation}
+With $a_i$ is the "Modal A" parameter linked to mode i.
+
+#+begin_src matlab
+  modal_a_M = diag(complex(modal_a(:, 1), modal_a(:, 2)));
+  modal_b_M = diag(complex(modal_b(:, 1), modal_b(:, 2)));
+
+  modal_a_ext_M = blkdiag(modal_a_M, conj(modal_a_M));
+  modal_b_ext_M = blkdiag(modal_b_M, conj(modal_b_M));
+#+end_src
 
 * Obtained Mode Shapes animations
   <<sec:mode_shape_display>>
@@ -441,6 +501,133 @@ 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.
 
+* Verify the validity of the Modal Model
+  <<sec:validity_measurement>>
+There are two main ways to verify the validity of the modal model
+- Synthesize FRF measurements that has been used to generate the modal model and compare
+- Synthesize FRF that has not yet been measured. Then measure that FRF and compare
+
+** Theory
+From the modal model, we want to synthesize the Frequency Response Functions that has been used to build the modal model.
+
+Let's recall that:
+- $M$ is the number of measured DOFs ($23 \times 3 = 69$)
+- $N$ is the number of modes identified ($21$)
+
+We then have that the FRF matrix $[H_{\text{syn}}]$ can be synthesize using the following formula:
+#+begin_important
+\begin{equation}
+  [H_{\text{syn}}(\omega)]_{M\times M} = [\Phi]_{M\times2N} \left[\frac{Q_r}{j\omega - s_r}\right]_{2N\times2N} [\Phi]_{2N\times M}^T
+\end{equation}
+with $Q_r = 1/M_{A_r}$
+#+end_important
+
+An alternative formulation is:
+\[ H_{pq}(s_i) = \sum_{r=1}^N \frac{A_{pqr}}{s_i - \lambda_r} + \frac{A_{pqr}^*}{s_i - \lambda_r^*} \]
+with:
+- $A_{pqr} = \frac{\psi_{pr}\psi_{qr}}{M_{A_r}}$, $M_{A_r}$ is called "Modal A"
+- $\psi_{pr}$: scaled modal coefficient for output DOF $p$, mode $r$
+- $\lambda_r$: complex modal frequency
+
+** Matlab Implementation
+#+begin_src matlab
+  Hsyn = zeros(69, 69, 801);
+
+  for i = 1:length(freqs)
+    Hsyn(:, :, i) = eigen_vec_ext_M*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M.';
+  end
+#+end_src
+
+Because the synthesize frequency response functions are representing the displacement response in $[m/N]$, we multiply each element of the FRF matrix by $(j \omega)^2$ in order to obtain the acceleration response in $[m/s^2/N]$.
+
+#+begin_src matlab
+  for i = 1:size(Hsyn, 1)
+    Hsyn(i, :, :) = squeeze(-Hsyn(i, :, :)).*(j*2*pi*freqs).^2;
+  end
+#+end_src
+
+** Original and Synthesize FRF matrix comparison
+#+begin_src matlab :exports none
+  acc_o = 15; dir_o = 1; dir_i = 1;
+
+  figure;
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(freqs, abs(squeeze(FRFs(3*(acc_o-1)+dir_o,          dir_i, :))), 'DisplayName', 'Original');
+  plot(freqs, abs(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :))), 'DisplayName', 'Synthesize');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+  legend('location', 'northwest');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(freqs, mod(180/pi*phase(squeeze(FRFs(3*(acc_o-1)+dir_o,          dir_i, :)))+180, 360)-180);
+  plot(freqs, mod(180/pi*phase(squeeze(Hsyn(3*(acc_o-1)+dir_o, 3*(11-1)+dir_i, :)))+180, 360)-180);
+  hold off;
+  yticks(-360:90:360); ylim([-180, 180]);
+  set(gca, 'xscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/compare_synthesize_original_frf.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:compare_synthesize_original_frf
+#+CAPTION: Comparison of the Original and Synthesize FRF matrix element
+[[file:figs/compare_synthesize_original_frf.png]]
+
+** Synthesize FRF that has not yet been measured
+#+begin_src matlab :exports none
+  dir_names = {'X', 'Y', 'Z'};
+
+  accs = [1]; dirs = [1:3];
+
+  figure;
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for acc_i = accs
+    for dir_i = dirs
+      plot(freqs, abs((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :))), 'DisplayName', sprintf('Acc %i - %s', acc_i, dir_names{dir_i}));
+    end
+  end
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude [$\frac{m}{N}$]');
+  legend('location', 'southwest');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for acc_i = accs
+    for dir_i = dirs
+      plot(freqs, mod(180/pi*phase((1./(j*2*pi*freqs').^2).*squeeze(Hsyn(3*(acc_i-1)+dir_i, 3*(acc_i-1)+dir_i, :)))+180, 360)-180);
+    end
+  end
+  hold off;
+  yticks(-360:90:360); ylim([-180, 180]);
+  set(gca, 'xscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/synthesize_frf_new_meas.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:synthesize_frf_new_meas
+#+CAPTION: Synthesize FRF that has not been measured
+[[file:figs/synthesize_frf_new_meas.png]]
+
 * Modal Complexity
   <<sec:modal_complexity>>
 A method of displaying *modal complexity* is by plotting the elements of the eigenvector on an *Argand diagram* (complex plane), such as the ones shown in figure [[fig:modal_complexity_small]].
@@ -527,10 +714,59 @@ The complexity of all the modes are compared on figure [[fig:modal_complexities]
 #+CAPTION: Modal complexity for each mode
 [[file:figs/modal_complexities.png]]
 
-* BKMK From local coordinates to global coordinates for the mode shapes
+* Modal Assurance Criterion (MAC)
+The MAC is calculated as the normalized scalar product of the two sets of vectors $\{\psi_A\}$ and $\{\psi_X\}$.
+The resulting scalars are arranged into the MAC matrix cite:pastor12_modal_assur_criter:
+\begin{equation}
+  \text{MAC}(r, q) = \frac{\left| \{\psi_A\}_r^T \{\psi_X\}_q^* \right|^2}{\left( \{\psi_A\}_r^T \{\psi_A\}_r^* \right) \left( \{\psi_X\}_q^T \{\psi_X\}_q^* \right)}
+\end{equation}
+
+An equivalent formulation is:
+\begin{equation}
+  \text{MAC}(r, q) = \frac{\left| \sum_{j=1}^n \{\psi_A\}_j \{\psi_X\}_j^* \right|^2}{\left( \sum_{j=1}^n |\{\psi_A\}_j|^2 \right) \left( \sum_{j=1}^n |\{\psi_X\}_j|^2 \right)}
+\end{equation}
+
+The MAC takes value between 0 (representing no consistent correspondence) and 1 (representing a consistent correspondence).
+
+We compute the autoMAC matrix that compares all the possible combinations of mode shape pairs for only one set of mode shapes. The result is shown on figure [[fig:automac]].
+
+#+begin_src matlab
+  autoMAC = eye(size(eigen_vec_M, 2));
+
+  for r = 1:size(eigen_vec_M, 2)
+    for q = r+1:size(eigen_vec_M, 2)
+      autoMAC(r, q) = abs(eigen_vec_M(r, :)*eigen_vec_M(q, :)')^2/((eigen_vec_M(r, :)*eigen_vec_M(r, :)')*(eigen_vec_M(q, :)*eigen_vec_M(q, :)'));
+      autoMAC(q, r) = autoMAC(r, q);
+    end
+  end
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  imagesc(autoMAC);
+  colormap('parula');
+  colorbar;
+  xticks(1:1:size(eigen_val_M, 1)); yticks(1:1:size(eigen_val_M, 1));
+  xticklabels(num2str(round(imag(diag(eigen_val_M))/2/pi, 1)));
+  xtickangle(90);
+  yticklabels(num2str(round(imag(diag(eigen_val_M))/2/pi, 1)));
+  set(gca,'YDir','normal');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/automac.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:automac
+#+CAPTION: AutoMAC
+[[file:figs/automac.png]]
+
+* From accelerometer coordinates to global coordinates for the mode shapes
   <<sec:global_mode_shapes>>
+
 ** Mathematical description
-#+begin_src latex :file local_to_global_coordinates.pdf :post pdf2svg(file=*this*, ext="png") :exports results
+#+begin_src latex :file local_to_global_coordinates_modes.pdf :post pdf2svg(file=*this*, ext="png") :exports results
   \newcommand\irregularcircle[2]{% radius, irregularity
     \pgfextra {\pgfmathsetmacro\len{(#1)+rand*(#2)}}
     +(0:\len pt)
@@ -540,7 +776,7 @@ The complexity of all the modes are compared on figure [[fig:modal_complexities]
     } -- cycle
   }
   \begin{tikzpicture}
-    \draw[rounded corners=1mm] (0, 0) \irregularcircle{3cm}{1mm};
+    \draw[rounded corners=1mm, fill=blue!30!white] (0, 0) \irregularcircle{3cm}{1mm};
 
     \node[] (origin) at (4, -1) {$\bullet$};
     \begin{scope}[shift={(origin)}]
@@ -557,22 +793,24 @@ The complexity of all the modes are compared on figure [[fig:modal_complexities]
     \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$};
+    \draw[->] (p1)node[]{$\bullet$}node[above]{$p_1$} -- ++( 1, 0.5)node[right]{$\xi_1$};
+    \draw[->] (p2)node[]{$\bullet$}node[above]{$p_2$} -- ++(-0.5, 1)node[right]{$\xi_2$};
+    \draw[->] (p3)node[]{$\bullet$}node[above]{$p_3$} -- ++( 1, 0.5)node[right]{$\xi_3$};
+    \draw[->] (p4)node[]{$\bullet$}node[above]{$p_4$} -- ++( 0.5, 1)node[right]{$\xi_4$};
   \end{tikzpicture}
 #+end_src
 
 #+RESULTS:
-[[file:figs/local_to_global_coordinates.png]]
+[[file:figs/local_to_global_coordinates_modes.png]]
+
+If we note $\xi_i = \{ \xi_{i,x}\ \xi_{i,y}\ \xi_{i,z} \}^T$ the mode shape corresponding to the accelerometer $i$.
 
 From the figure above, 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
+  \xi_1 &= \xi + \Omega p_1\\
+  \xi_2 &= \xi + \Omega p_2\\
+  \xi_3 &= \xi + \Omega p_3\\
+  \xi_4 &= \xi + \Omega p_4
 \end{align*}
 
 With
@@ -583,30 +821,13 @@ With
   -\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:
+$\xi$ and $\Omega$ represent to mode shapes in translation and rotation of the solid expressed in the frame $\{O\}$.
 
-\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\}$:
+We can rearrange the equations in a matrix form, and then we obtain the mode shape 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
+  \xi_x \\ \xi_y \\ \xi_z \\ \hline \Omega_x \\ \Omega_y \\ \Omega_z
 \end{bmatrix} =
 \left[\begin{array}{ccc|ccc}
   1 & 0 & 0 & 0 & p_{1z} & -p_{1y} \\
@@ -617,74 +838,181 @@ and then we obtain the velocity and rotation of the solid in the wanted frame $\
   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}
+  \xi_{1,x} \\ \xi_{1,y} \\ \xi_{1,z} \\\hline \vdots \\\hline \xi_{4,x} \\ \xi_{4,y} \\ \xi_{4,z}
 \end{bmatrix}
 \end{equation}
 
 This inversion is equivalent to a mean square problem.
 
 ** Matlab Implementation
+The obtained mode shapes matrix that gives the mode shapes of each solid bodies with respect to the fixed frame $\{O\}$, =mode_shapes_O=, is an $n \times p \times q$ with:
+- $n$ is the number of modes: 21
+- $p$ is the number of DOFs for each solid body: 6
+- $q$ is the number of solid bodies: 6
+
 #+begin_src matlab
-  mode_shapes_O = zeros(mod_n, length(solid_names), 6);
+  mode_shapes_O = zeros(mod_n, 6, length(solid_names));
 
   for mod_i = 1:mod_n
     for solid_i = 1:length(solid_names)
-      solids_i = solids.(solid_names{solid_i});
+      acc_i = solids.(solid_names{solid_i});
 
-      Y = reshape(squeeze(modes(mod_i, solids_i, :))', [], 1);
+      Y = mode_shapes(mod_i, :, acc_i);
+      Y = Y(:);
 
-      A = zeros(3*length(solids_i), 6);
-      for i = 1:length(solids_i)
+      A = zeros(3*length(acc_i), 6);
+      for i = 1:length(acc_i)
         A(3*(i-1)+1:3*i, 1:3) = eye(3);
 
-        A(3*(i-1)+1:3*i, 4:6) = [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];
+        A(3*(i-1)+1:3*i, 4:6) = [ 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];
       end
 
-      mode_shapes_O(mod_i, solid_i, :) = A\Y;
+      mode_shapes_O(mod_i, :, solid_i) = A\Y;
     end
   end
 #+end_src
 
-** Modal Matrices
-We arrange the obtained modal parameters into matrices:
-\[ \Omega = \begin{bmatrix}
-  \omega_1^2 &        & 0 \\
-             & \ddots &   \\
-  0          &        & \omega_n^2
-\end{bmatrix}; \quad \Psi = \begin{bmatrix}
+We then rearrange the eigen vectors in another way:
+\[ \Psi_O = \begin{bmatrix}
   & & \\
-  \{\psi_1\} & \dots & \{\psi_n\} \\
+  \{\psi_1\} & \dots & \{\psi_N\} \\
   & &
-\end{bmatrix} \]
+\end{bmatrix}_{M \times N} \]
 with
-\[ \{\psi_1\} = \begin{Bmatrix} \psi_{1_x} & \psi_{2_x} & \dots & \psi_{6_x} & \psi_{1_y} & \dots & \psi_{1\Omega_x} & \dots & \psi_{6\Omega_z} \end{Bmatrix}^T \]
+\[ \{\psi\}_r = \begin{Bmatrix}
+  \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z}
+\end{Bmatrix}^T \]
+
+With $M = 6 \times 6$ is the new number of DOFs and $N=21$ is the number of modes.
+
+#+begin_src matlab
+  eigen_vec_O = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)]).';
+
+  eigen_vec_ext_O = [eigen_vec_O, conj(eigen_vec_O)];
+#+end_src
+
+** FRF matrix synthesis
+*** Mathematical description
+The synthesize FRF matrix =Hsyn_O= is an $n \times n \times q$ with:
+- $n$ is the number of DOFs of the considered 6 solid-bodies: $6 \times 6 = 36$
+- $q$ is the number of frequency points $\omega_i$
+
+For each frequency point $\omega_i$, the FRF matrix =Hsyn_O= is a $n\times n$ matrix:
+\begin{equation}
+  \text{H}_{\text{syn, O}}(\omega_i) = \begin{bmatrix}
+  \frac{D_{1,T_x}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,T_x}}{F_{1, T_y}}(\omega_i) & \dots & \frac{D_{1,T_x}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{1,T_y}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,T_y}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{1,T_y}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{1,T_z}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,T_z}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{1,T_z}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{1,R_x}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,R_x}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{1,R_x}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{1,R_y}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,R_y}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{1,R_y}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{1,R_z}}{F_{1, T_x}}(\omega_i) & \frac{D_{1,R_z}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{1,R_z}}{F_{6, R_z}}(\omega_i) \\
+  \frac{D_{2,T_x}}{F_{1, T_x}}(\omega_i) & \frac{D_{2,T_x}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{2,T_x}}{F_{6, R_z}}(\omega_i) \\
+  \vdots                                 & \vdots                                 &       & \vdots                                 \\
+  \frac{D_{6,R_z}}{F_{1, T_x}}(\omega_i) & \frac{D_{6,R_z}}{F_{1, T_y}}(\omega_i) &       & \frac{D_{6,R_z}}{F_{6, R_z}}(\omega_i)
+\end{bmatrix}
+\end{equation}
+where $D_i$ corresponds to the solid body number i, $F_i$ is a force/torque applied on the solid body number i at the center of the frame $\{O\}$.
 
 #+begin_warning
-  How to add damping to the eigen matrices?
+  As the mode shapes are expressed with respect to the frame $\{O\}$, the obtain FRF matrix are the responses of the system from forces and torques applied to the solid bodies at the origin of $\{O\}$ to the displacement of the solid bodies with respect to the frame $\{O\}$.
 #+end_warning
 
+*** Matlab Implementation
 #+begin_src matlab
-  eigen_value_M = diag(freqs*2*pi);
-  eigen_vector_M = reshape(mode_shapes_O, [mod_n, 6*length(solid_names)])';
+  Hsyn_O = zeros(36, 36, 801);
+
+  for i = 1:length(freqs)
+    Hsyn_O(:, :, i) = eigen_vec_ext_O*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_O.';
+  end
 #+end_src
 
-** TODO Normalization of mode shapes?
-We normalize each column of the eigen vector matrix.
-
-Then, each eigenvector as a norm of 1.
+Because the synthesize frequency response functions are representing the displacement response in $[m/N]$, we multiply each element of the FRF matrix by $(j \omega)^2$ in order to obtain the acceleration response in $[m/s^2/N]$.
 
 #+begin_src matlab
-  eigen_vector_M = eigen_vector_M./vecnorm(eigen_vector_M);
+  for i = 1:size(Hsyn_O, 1)
+    Hsyn_O(i, :, :) = -squeeze(Hsyn_O(i, :, :)).*(j*2*pi*freqs).^2;
+  end
 #+end_src
 
-#+begin_important
-  Should we do such normalization?
-#+end_important
+*** TODO Test to have the original inputs
+#+begin_src matlab
+  Hsyn_O = zeros(36, 3, 801);
+
+  for i = 1:length(freqs)
+    Hsyn_O(:, :, i) = eigen_vec_ext_O*(inv(modal_a_ext_M)/(diag(j*2*pi*freqs(i) - diag(eigen_val_ext_M))))*eigen_vec_ext_M(10*3+1:11*3, :).';
+  end
+
+  for i = 1:size(Hsyn_O, 1)
+    Hsyn_O(i, :, :) = -squeeze(Hsyn_O(i, :, :)).*(j*2*pi*freqs).^2;
+  end
+#+end_src
+
+** TODO Verification
+#+begin_src matlab :exports none
+  solid_o = 6; dir_o = 3;
+  solid_i = 3; dir_i = 3;
+
+  figure;
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(freqs, abs(squeeze(FRFs_O((solid_o-1)*6+dir_o, dir_i, :))), 'DisplayName', 'Original');
+  plot(freqs, abs(squeeze(Hsyn_O((solid_o-1)*6+dir_o, dir_i, :))), 'DisplayName', 'Synthesize');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+  legend('location', 'northwest');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(freqs, mod(180/pi*phase(squeeze(FRFs_O((solid_o-1)*6+dir_o, dir_i, :)))+180, 360)-180);
+  plot(freqs, mod(180/pi*phase(squeeze(Hsyn_O((solid_o-1)*6+dir_o, dir_i, :)))+180, 360)-180);
+  hold off;
+  yticks(-360:90:360);
+  set(gca, 'xscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+** New FRFs
+#+begin_src matlab :exports none
+  dof_names = {'X', 'Y', 'Z', '\theta_X', '\theta_Y', '\theta_Z'};
+  solid_i = 6; dirs = [1:3];
+
+  figure;
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for dir_i = dirs
+    plot(freqs, abs(squeeze(Hsyn_O((solid_i-1)*6+dir_i, (solid_i-1)*6+dir_i, :))), 'DisplayName', sprintf('Solid %s - %s', solid_names{solid_i}, dof_names{dir_i}));
+  end
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude [$\frac{m/s^2}{N}$]');
+  legend('location', 'northwest');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for dir_i = dirs
+    plot(freqs, mod(180/pi*phase(squeeze(Hsyn_O((solid_i-1)*6+dir_i, (solid_i-1)*6+dir_i, :)))+180, 360)-180);
+  end
+  hold off;
+  yticks(-360:90:360);
+  set(gca, 'xscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
 
 * 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.
+They should have the same size $n \times n$ where $n$ is the number of mode_shapes 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.
@@ -759,4 +1087,232 @@ If we look at the z motion for instance, we will find that we cannot neglect tha
   test = mode_shapes_O(10, 2, :)/norm(squeeze(mode_shapes_O(10, 2, :)));
 #+end_src
 
-* TODO Synthesis of FRF curves from the modal parameters
+* TODO Mode shapes expressed in a frame at the CoM of the solid bodies
+#+begin_warning
+  Maybe the mass-normalized eigen vectors can only be obtained for undamped systems (real eigen values and vectors).
+#+end_warning
+
+** Modes shapes transformation
+First, we load the position of the Center of Mass of each solid body with respect to the point of interest.
+#+begin_src matlab
+  model_com = table2array(readtable('mat/model_solidworks_com.txt', 'ReadVariableNames', false));
+#+end_src
+
+Then, we reshape the data.
+#+begin_src matlab
+  model_com = reshape(model_com , [3, 6]);
+#+end_src
+
+Now we convert the mode shapes expressed in the DOF of the accelerometers to the DoF of each solid body (translations and rotations) with respect to their own CoM.
+#+begin_src matlab
+  mode_shapes_CoM = zeros(mod_n, 6, length(solid_names));
+
+  for mod_i = 1:mod_n
+    for solid_i = 1:length(solid_names)
+      acc_i = solids.(solid_names{solid_i});
+
+      Y = mode_shapes(mod_i, :, acc_i);
+      Y = Y(:);
+
+      A = zeros(3*length(acc_i), 6);
+      for i = 1:length(acc_i)
+        acc_pos_com = acc_pos(i, :).' - model_com(:, i);
+
+        A(3*(i-1)+1:3*i, 1:3) = eye(3);
+        A(3*(i-1)+1:3*i, 4:6) = [ 0               acc_pos_com(3) -acc_pos_com(2) ;
+                                 -acc_pos_com(3)  0               acc_pos_com(1) ;
+                                  acc_pos_com(2) -acc_pos_com(1)  0];
+      end
+
+      mode_shapes_CoM(mod_i, :, solid_i) = A\Y;
+    end
+  end
+#+end_src
+
+We then rearrange the eigen vectors in another way:
+\[ \Psi_{\text{CoM}} = \begin{bmatrix}
+  & & \\
+  \{\psi_1\} & \dots & \{\psi_N\} \\
+  & &
+\end{bmatrix}_{M \times N} \]
+with
+\[ \{\psi\}_r = \begin{Bmatrix}
+  \psi_{1, x} & \psi_{1, y} & \psi_{1, z} & \psi_{1, \theta_x} & \psi_{1, \theta_y} & \psi_{1, \theta_z} & \psi_{2, x} & \dots & \psi_{6, \theta_z}
+\end{Bmatrix}^T \]
+
+With $M = 6 \times 6$ is the new number of DOFs and $N=21$ is the number of modes.
+
+Each eigen vector is normalized.
+
+#+begin_src matlab
+  eigen_vec_CoM = reshape(mode_shapes_CoM, [mod_n, 6*length(solid_names)]).';
+  % eigen_vec_CoM = eigen_vec_CoM./vecnorm(eigen_vec_CoM);
+  eigen_vec_ext_CoM = [eigen_vec_CoM, conj(eigen_vec_CoM)];
+#+end_src
+
+** Mass matrix from solidworks model
+#+begin_src bash :results none
+  sed '/Mass = /!d' mat/1_granite_bot.txt | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >  mat/model_solidworks_mass.txt
+  sed '/Mass = /!d' mat/2_granite_top.txt | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
+  sed '/Mass = /!d' mat/3_ty.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
+  sed '/Mass = /!d' mat/4_ry.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
+  sed '/Mass = /!d' mat/5_rz.txt          | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
+  sed '/Mass = /!d' mat/6_hexa.txt        | sed 's/^\s*\t*Mass = //' | sed 's/\([0-9.]*\) kilograms\r/\1/' >> mat/model_solidworks_mass.txt
+#+end_src
+
+
+#+begin_src bash :results none
+  sed '/I[xyz][xyz]/!d' mat/1_granite_bot.txt | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >  mat/model_solidworks_inertia.txt
+  sed '/I[xyz][xyz]/!d' mat/2_granite_top.txt | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
+  sed '/I[xyz][xyz]/!d' mat/3_ty.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
+  sed '/I[xyz][xyz]/!d' mat/4_ry.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
+  sed '/I[xyz][xyz]/!d' mat/5_rz.txt          | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
+  sed '/I[xyz][xyz]/!d' mat/6_hexa.txt        | sed 's/^\s*\t*//' | sed 's/I[xyz][xyz] = //g' | sed $'s/\r//' >> mat/model_solidworks_inertia.txt
+#+end_src
+
+#+begin_src matlab
+  model_mass = table2array(readtable('mat/model_solidworks_mass.txt', 'ReadVariableNames', false));
+  model_inertia = table2array(readtable('mat/model_solidworks_inertia.txt', 'ReadVariableNames', false));
+#+end_src
+
+\[ M = \begin{bmatrix}
+  M_1 &     &     &     &     &     \\
+      & M_2 &     &     &     &     \\
+      &     & M_3 &     &     &     \\
+      &     &     & M_4 &     &     \\
+      &     &     &     & M_5 &     \\
+      &     &     &     &     & M_6 \\
+\end{bmatrix}, \text{ with } M_i = \begin{bmatrix}
+  m_i &     &     &          &          &          \\
+      & m_i &     &          &          &          \\
+      &     & m_i &          &          &          \\
+      &     &     & I_{i,xx} & I_{i,yx} & I_{i,zx} \\
+      &     &     & I_{i,xy} & I_{i,yy} & I_{i,zy} \\
+      &     &     & I_{i,xz} & I_{i,yz} & I_{i,zz} \\
+\end{bmatrix} \]
+
+#+begin_src matlab
+  M = zeros(6*6, 6*6);
+
+  for i = 1:6
+    M((i-1)*6+1:i*6, (i-1)*6+1:i*6) = blkdiag(model_mass(i)*eye(3), model_inertia((i-1)*3+1:i*3, 1:3));
+  end
+#+end_src
+
+** Mass-normalized Eigen Vectors
+To do so, it seems that we need to know the mass matrix $[M]$. Then:
+\[ \{\phi\}_r = \frac{1}{\sqrt{m_r}} \{\psi\}_r \text{ where } m_r = \{\psi\}_r^T [M] \{\psi\}_r \]
+
+Mass-normalized eigen vectors are very important for the synthesis and spatial model extraction.
+
+Mass Matrix could be estimated from the SolidWorks model.
+
+We compute the modal masses that will be used for normalization.
+#+begin_src matlab
+  mr = zeros(size(eigen_vector_matrix_CoM, 2), 1);
+
+  for i = 1:length(mr)
+    mr(i) = real(eigen_vector_matrix_CoM(:, i).' * M * eigen_vector_matrix_CoM(:, i));
+  end
+#+end_src
+
+#+begin_src matlab
+  figure;
+  plot(freqs_m, mr, 'ko');
+  xlabel('Mode Frequency [Hz]'); ylabel('Modal mass [kg]');
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/modal_mass.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:modal_mass
+#+CAPTION: Modal masses
+[[file:figs/modal_mass.png]]
+
+And finally, we compute the mass-normalized eigen vectors.
+#+begin_src matlab
+  eigen_vector_mass_CoM = zeros(size(eigen_vector_matrix_CoM));
+
+  for i = 1:size(eigen_vector_matrix_CoM, 2)
+    eigen_vector_mass_CoM(:, i) = 1/sqrt(mr(i)) * eigen_vector_matrix_CoM(:, i);
+  end
+#+end_src
+
+Verification that
+\[ [\Phi]^T [M] [\Phi] = [I] \]
+
+#+begin_src matlab
+  eigen_vector_matrix_CoM.'*M*eigen_vector_matrix_CoM
+#+end_src
+
+#+begin_src matlab
+  eigen_vector_mass_CoM*M*eigen_vector_mass_CoM.'
+#+end_src
+
+Other test for normalized eigen vectors
+#+begin_src matlab
+  eigen_vector_mass_CoM = (eigen_vector_matrix_CoM.'*diag(diag(M))*eigen_vector_matrix_CoM)^(-0.5) * eigen_vector_matrix_CoM';
+  eigen_vector_mass_CoM = eigen_vector_mass_CoM.';
+#+end_src
+
+** Full Response Model from modal model (synthesis)
+In general, the form of response model with which we are concerned is an *FRF matrix* whose order is dictated by the number of coordinates $n$ included in the test.
+Also, as explained, it is normal in practice to measured and to analyze just a *subset of the FRF matrix* but rather to measure the full FRF matrix.
+Usually *one column* or *one row* with a few additional elements are measured.
+
+Thus, if we are to construct an acceptable response model, it will be necessary to synthesize those elements which have not been directly measured.
+However, in principle, this need present no major problem as it is possible to compute the full FRF matrix from a modal model using:
+\begin{equation}
+  [H]_{n\times n} = [\Phi]_{n\times m} [\lambda_r^2 - \omega^2]_{m\times m}^{-1} [\Phi]_{m\times n}^T
+\end{equation}
+
+$\{\Phi\}$ is a *mass-normalized* eigen vector.
+
+#+begin_src matlab
+  FRF_matrix_CoM = zeros(size(eigen_vector_mass_CoM, 1), size(eigen_vector_mass_CoM, 1), length(freqs));
+
+  for i = 1:length(freqs)
+    FRF_matrix_CoM(:, :, i) = eigen_vector_mass_CoM*inv(eigen_value_M^2 - (2*pi*freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass_CoM.';
+  end
+#+end_src
+
+#+begin_src matlab
+  exc_dir = 3;
+  meas_mass = 6;
+  meas_dir = 3;
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(FRF_matrix_CoM((meas_mass-1)*6 + meas_dir, 6*2+exc_dir, :))));
+  plot(freqs, abs(squeeze(FRFs_CoM((meas_mass-1)*6 + meas_dir, exc_dir, :))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude');
+  xlim([1, 200]);
+#+end_src
+
+#+begin_src matlab
+  FRF_matrix = zeros(size(eigen_vector_mass, 1), size(eigen_vector_mass, 1), length(freqs));
+
+  for i = 1:length(freqs)
+    FRF_matrix(:, :, i) = eigen_vector_mass*inv(eigen_value_M - (freqs(i)*eye(size(eigen_value_M, 1)))^2)*eigen_vector_mass.';
+  end
+#+end_src
+
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(FRFs_O(1, 1, :))));
+  plot(freqs, abs(squeeze(FRF_matrix(1, 1, :))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude');
+  xlim([1, 200]);
+#+end_src
+
+* TODO Residues
+* Bibliography                                                       :ignore:
+bibliographystyle:unsrt
+bibliography:ref.bib

modal-analysis/ref.bib

@@ -0,0 +1,38 @@
+@article{wang11_extrac_real_modes_physic_matric,
+  author =       {Tong Wang and Lingmi Zhang and Kong Fah Tee},
+  title =        {Extraction of Real Modes and Physical Matrices From Modal
+                  Testing},
+  journal =      {Earthquake Engineering and Engineering Vibration},
+  volume =       10,
+  number =       2,
+  pages =        {219-227},
+  year =         2011,
+  doi =          {10.1007/s11803-011-0060-6},
+  url =          {https://doi.org/10.1007/s11803-011-0060-6},
+  DATE_ADDED =   {Tue Jul 9 15:51:21 2019},
+}
+
+@article{pastor12_modal_assur_criter,
+  author = {Miroslav Pastor and Michal Binda and Tom{\'a}{\v{s}} Har{\v{c}}arik},
+  title = {Modal Assurance Criterion},
+  journal = {Procedia Engineering},
+  volume = {48},
+  number = {nil},
+  pages = {543-548},
+  year = {2012},
+  doi = {10.1016/j.proeng.2012.09.551},
+  url = {https://doi.org/10.1016/j.proeng.2012.09.551},
+  DATE_ADDED = {Thu Jul 11 13:51:57 2019},
+}
+
+@book{ewins00_modal,
+  author = {Ewins, DJ},
+  title = {Modal testing: theory, practice and application},
+  year = {2000},
+  publisher = {Wiley-Blackwell},
+  journal = {Research studies Pre, 2nd ed., ISBN-13},
+  pages = {978--0863802188},
+  isbn = {0863802184},
+  address = {Baldock, Hertfordshire, England Philadelphia, PA},
+  keywords = {favorite, identification},
+}