Add importation of FRF and COH

modal-analysis/modes_analysis.org

@@ -347,7 +347,7 @@ The obtained area of this polygon is then compared with the area of the circle w
 A little complex mode 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]].
 
-#+begin_src matlab :export none
+#+begin_src matlab :exports none
   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)));
@@ -373,7 +373,7 @@ The complexity of all the modes are compared on figure [[fig:modal_complexities]
 #+CAPTION: Modal Complexity of one mode with small complexity
 [[file:figs/modal_complexity_small.png]]
 
-#+begin_src matlab :export none
+#+begin_src matlab :exports none
   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)));
@@ -399,7 +399,7 @@ The complexity of all the modes are compared on figure [[fig:modal_complexities]
 #+CAPTION: Modal Complexity of one higly complex mode
 [[file:figs/modal_complexity_high.png]]
 
-#+begin_src matlab :export none
+#+begin_src matlab :exports none
   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)));
@@ -508,4 +508,261 @@ 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
 
+* Importation of measured FRF curves
+There are 24 measurements files corresponding to 24 series of impacts:
+- 3 directions, 8 sets of 3 accelerometers
+For each measurement file, the FRF and coherence between the impact and the 9 accelerations measured.
+
+In reality: 4 sets of 10 things
+
+#+begin_src matlab
+  a = load('./modal_analysis/frf_coh/Measurement1.mat');
+#+end_src
+
+#+begin_src matlab
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Mod)
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+  title(sprintf('From %s, to %s', FFT1_AvXSpc_2_1_RfName, FFT1_AvXSpc_2_1_RpName))
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(a.FFT1_AvXSpc_2_1_RMS_X_Val, a.FFT1_AvXSpc_2_1_RMS_Y_Phas)
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+* Importation of measured FRF curves to global FRF matrix
+FRF matrix $n \times p$:
+- $n$ is the number of measurements: $3 \times 24$
+- $p$ is the number of excitation inputs: 3
+
+23 measurements: 3 accelerometers
+
+\begin{equation}
+  \text{FRF}(\omega_i) = \begin{bmatrix}
+  \frac{D_{1_x}}{F_x}(\omega_i) & \frac{D_{1_x}}{F_y}(\omega_i) & \frac{D_{1_x}}{F_z}(\omega_i) \\
+  \frac{D_{1_y}}{F_x}(\omega_i) & \frac{D_{1_y}}{F_y}(\omega_i) & \frac{D_{1_y}}{F_z}(\omega_i) \\
+  \frac{D_{1_z}}{F_x}(\omega_i) & \frac{D_{1_z}}{F_y}(\omega_i) & \frac{D_{1_z}}{F_z}(\omega_i) \\
+  \frac{D_{2_x}}{F_x}(\omega_i) & \frac{D_{2_x}}{F_y}(\omega_i) & \frac{D_{2_x}}{F_z}(\omega_i) \\
+  \vdots              & \vdots              & \vdots              \\
+  \frac{D_{23_z}}{F_x}(\omega_i) & \frac{D_{23_z}}{F_y}(\omega_i) & \frac{D_{23_z}}{F_z}(\omega_i) \\
+\end{bmatrix}
+\end{equation}
+
+#+begin_src matlab
+  n_meas = 24;
+  n_acc = 23;
+
+  dirs = 'XYZ';
+
+  % Number of Accelerometer * DOF for each acccelerometer / Number of excitation / frequency points
+  FRFs = zeros(3*n_acc, 3, 801);
+  COHs = zeros(3*n_acc, 3, 801);
+
+  % Loop through measurements
+  for i = 1:n_meas
+    % Load the measurement file
+    meas = load(sprintf('./modal_analysis/frf_coh/Measurement%i.mat', i));
+
+    % First: determine what is the exitation (direction and sign)
+    exc_dir = meas.FFT1_AvXSpc_2_1_RMS_RfName(end);
+    exc_sign = meas.FFT1_AvXSpc_2_1_RMS_RfName(end-1);
+    % Determine what is the correct excitation sign
+    exc_factor = str2num([exc_sign, '1']);
+    if exc_dir ~= 'Z'
+      exc_factor = exc_factor*(-1);
+    end
+
+    % Then: loop through the nine measurements and store them at the correct location
+    for j = 2:10
+      % Determine what is the accelerometer and direction
+      [indices_acc_i] = strfind(meas.(sprintf('FFT1_H1_%i_1_RpName', j)), '.');
+      acc_i = str2num(meas.(sprintf('FFT1_H1_%i_1_RpName', j))(indices_acc_i(1)+1:indices_acc_i(2)-1));
+
+      meas_dir  = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end);
+      meas_sign = meas.(sprintf('FFT1_H1_%i_1_RpName', j))(end-1);
+      % Determine what is the correct measurement sign
+      meas_factor = str2num([meas_sign, '1']);
+      if meas_dir ~= 'Z'
+        meas_factor = meas_factor*(-1);
+      end
+
+      % FRFs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
+      % COHs(acc_i+n_acc*(find(dirs==meas_dir)-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
+
+      FRFs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = exc_factor*meas_factor*meas.(sprintf('FFT1_H1_%i_1_Y_ReIm', j));
+      COHs(find(dirs==meas_dir)+3*(acc_i-1), find(dirs==exc_dir), :) = meas.(sprintf('FFT1_Coh_%i_1_RMS_Y_Val', j));
+    end
+  end
+  freqs = meas.FFT1_Coh_10_1_RMS_X_Val;
+#+end_src
+
+* Analysis of some FRFs
+#+begin_src matlab
+  acc_i = 3;
+  acc_dir = 1;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(freqs, abs(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(acc_i-1), exc_dir, :))), 360)-180);
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+
+#+begin_src matlab
+  figure;
+  hold on;
+  for i = 1:3*n_acc
+    plot(freqs, squeeze(COHs(i, 1, :)), 'color', [0, 0, 0, 0.2]);
+  end
+  hold off;
+  xlabel('Frequency [Hz]');
+  ylabel('Coherence [\%]');
+#+end_src
+
+
+Composite Response Function.
+
+We here sum the norm instead of the complex numbers.
+
+#+begin_src matlab
+  HHx = squeeze(sum(abs(FRFs(:, 1, :))));
+  HHy = squeeze(sum(abs(FRFs(:, 2, :))));
+  HHz = squeeze(sum(abs(FRFs(:, 3, :))));
+  HH  = squeeze(sum([HHx, HHy, HHz], 2));
+#+end_src
+
+#+begin_src matlab
+  exc_dir = 3;
+
+  figure;
+  hold on;
+  for i = 1:3*n_acc
+    plot(freqs, abs(squeeze(FRFs(i, exc_dir, :))), 'color', [0, 0, 0, 0.2]);
+  end
+  plot(freqs, abs(HHx));
+  plot(freqs, abs(HHy));
+  plot(freqs, abs(HHz));
+  plot(freqs, abs(HH), 'k');
+  hold off;
+  set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude');
+  xlim([1, 200]);
+#+end_src
+
+* From local coordinates to global coordinates with the FRFs
+#+begin_src matlab
+  % Number of Solids * DOF for each solid / Number of excitation / frequency points
+  FRFs_O = zeros(length(solid_names)*6, 3, 801);
+
+  for exc_dir = 1:3
+    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)
+        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
+
+      for i = 1:801
+        FRFs_O((solid_i-1)*6+1:solid_i*6, exc_dir, i) = A\FRFs((solids_i(1)-1)*3+1:solids_i(end)*3, exc_dir, i);
+      end
+    end
+  end
+#+end_src
+
+* Analysis of some FRF in the global coordinates
+#+begin_src matlab
+  solid_i = 6;
+  dir_i = 1;
+  exc_dir = 1;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+dir_i, exc_dir, :))), 360)-180);
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+* Compare global coordinates to local coordinates
+#+begin_src matlab
+  solid_i = 1;
+  acc_dir = 3;
+  exc_dir = 3;
+
+  figure;
+
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  for i = solids.(solid_names{solid_i})
+    plot(freqs, abs(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))));
+  end
+  plot(freqs, abs(squeeze(FRFs_O((solid_i-1)*6+acc_dir, exc_dir, :))), '-k');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Amplitude');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  for i = solids.(solid_names{solid_i})
+    plot(freqs, mod(180+180/pi*phase(squeeze(FRFs(acc_dir+3*(i-1), exc_dir, :))), 360)-180);
+  end
+  plot(freqs, mod(180+180/pi*phase(squeeze(FRFs_O((solid_i-1)*6+acc_dir, exc_dir, :))), 360)-180, '-k');
+  hold off;
+  ylim([-180, 180]); yticks(-180:90:180);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+  set(gca, 'xscale', 'log');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([1, 200]);
+#+end_src
+
+
 * TODO Synthesis of FRF curves