Add data analysis of marble dynamics

disturbance-sr-rz/index.org

@@ -30,6 +30,7 @@
 :END:
 
 * Experimental Setup
+*Setup*:
 All the stages are OFF.
 
 Two geophone are use:
@@ -43,9 +44,12 @@ Two voltage amplifiers are used, their setup is:
 
 A first order low pass filter is also added at the input of the voltage amplifiers.
 
+
 *Goal*:
 - Identify the vibrations induced by the rotation of the Slip-Ring and Spindle
 
+
+*Measurements*:
 Three measurements are done:
 | Measurement File   | Description                                                          |
 |--------------------+----------------------------------------------------------------------|

disturbance-ty-sr/index.org

@@ -1,4 +1,4 @@
-#+TITLE:Measurements
+#+TITLE:Vibrations induced by both the translation stage and the slip-ring
 :DRAWER:
 #+STARTUP: overview
 

dynamical-meas-granite/index.org

@@ -29,17 +29,281 @@
 #+PROPERTY: header-args:shell  :eval no-export
 :END:
 
+* Experimental Setup
+*Setup*:
+One geophone is located on the marble, the other one on the floor (see figure [[fig:experimental_setup]]).
 
-One geophone is located on the marble, the other one on the floor.
 
-col1: floor
+Two geophone are use:
-col2: marble
+- One on the floor (corresponding to the first column in the data)
+- One at the marble location (corresponding to the second column in the data)
 
-Each of the signal is amplified by voltage amplifier:
+Two voltage amplifiers are used, their setup is:
-- 60db
+- gain of 60dB
-- AC
+- AC/DC switch on AC
-- 1kHz
+- Low pass filter at 1kHz
 
-- meas 37: Z direction
+A first order low pass filter is also added at the input of the voltage amplifiers.
-- meas 38: N direction
+
-- meas 39: E direction
+*Goal*:
+- Identify the marble dynamics in all the directions
+
+
+*Measurements*:
+Three measurements are done:
+| Measurement File   | Description |
+|--------------------+-------------|
+| =mat/data_037.mat= | Z direction |
+| =mat/data_038.mat= | N direction |
+| =mat/data_039.mat= | E direction |
+
+Each of the measurement =mat= file contains one =data= array with 3 columns:
+| Column number | Description       |
+|---------------+-------------------|
+|             1 | Geophone - Floor  |
+|             2 | Geophone - Marble |
+|             3 | Time              |
+
+#+name: fig:experimental_setup
+#+caption: Picture of the experimental setup
+#+attr_html: :width 500px
+[[file:./img/IMG_20190513_161729.jpg]]
+
+#+name: fig:experimental_setup_bix
+#+caption: Picture of the experimental setup
+#+attr_html: :width 500px
+[[file:./img/IMG_20190513_161718.jpg]]
+
+* Data Analysis
+  :PROPERTIES:
+  :header-args:matlab+: :tangle matlab/marble_dynamics.m
+  :header-args:matlab+: :comments org :mkdirp yes
+  :END:
+  <<sec:marble_dynamics>>
+
+** ZIP file containing the data and matlab files                     :ignore:
+#+begin_src bash :exports none :results none
+  if [ matlab/marble_dynamics.m -nt data/marble_dynamics.zip ]; then
+    cp matlab/marble_dynamics.m marble_dynamics.m;
+    zip data/marble_dynamics \
+        mat/data_037.mat \
+        mat/data_038.mat \
+        mat/data_039.mat \
+        marble_dynamics.m
+    rm marble_dynamics.m;
+  fi
+#+end_src
+
+#+begin_note
+  All the files (data and Matlab scripts) are accessible [[file:data/marble_dynamics.zip][here]].
+#+end_note
+
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+<<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+<<matlab-init>>
+#+end_src
+
+** Load data
+#+begin_src matlab
+  m_z = load('mat/data_037.mat', 'data'); m_z = m_z.data;
+  m_n = load('mat/data_038.mat', 'data'); m_n = m_n.data;
+  m_e = load('mat/data_039.mat', 'data'); m_e = m_e.data;
+#+end_src
+
+** Time domain plots
+#+begin_src matlab
+  figure;
+  subplot(1, 3, 1);
+  hold on;
+  plot(m_z(:, 3), m_z(:, 2), 'DisplayName', 'Marble - Z');
+  plot(m_z(:, 3), m_z(:, 1), 'DisplayName', 'Floor - Z');
+  hold off;
+  xlabel('Time [s]'); ylabel('Voltage [V]');
+  xlim([0, 100]); ylim([-2 2]);
+  legend('Location', 'northeast');
+
+  subplot(1, 3, 2);
+  hold on;
+  plot(m_n(:, 3), m_n(:, 2), 'DisplayName', 'Marble - N');
+  plot(m_n(:, 3), m_n(:, 1), 'DisplayName', 'Floor - N');
+  hold off;
+  xlabel('Time [s]'); ylabel('Voltage [V]');
+  xlim([0, 100]); ylim([-2 2]);
+  legend('Location', 'northeast');
+
+  subplot(1, 3, 3);
+  hold on;
+  plot(m_e(:, 3), m_e(:, 2), 'DisplayName', 'Marble - E');
+  plot(m_e(:, 3), m_e(:, 1), 'DisplayName', 'Floor - E');
+  hold off;
+  xlabel('Time [s]'); ylabel('Voltage [V]');
+  xlim([0, 100]); ylim([-2 2]);
+  legend('Location', 'northeast');
+#+end_src
+
+#+NAME: fig:marble_floor_motion_time
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/marble_floor_motion_time.pdf" :var figsize="full-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:marble_floor_motion_time
+#+CAPTION: Floor and ground motion
+#+RESULTS: fig:marble_floor_motion_time
+[[file:figs/marble_floor_motion_time.png]]
+
+** Compute the power spectral densities
+We first compute some parameters that will be used for the PSD computation.
+#+begin_src matlab :results none
+  dt = m_z(2, 3)-m_z(1, 3);
+
+  Fs = 1/dt; % [Hz]
+
+  win = hanning(ceil(10*Fs));
+#+end_src
+
+Then we compute the Power Spectral Density using =pwelch= function.
+#+begin_src matlab
+  [px_fz, f] = pwelch(m_z(:, 1), win, [], [], Fs);
+  [px_gz, ~] = pwelch(m_z(:, 2), win, [], [], Fs);
+
+  [px_fn, ~] = pwelch(m_n(:, 1), win, [], [], Fs);
+  [px_gn, ~] = pwelch(m_n(:, 2), win, [], [], Fs);
+
+  [px_fe, ~] = pwelch(m_e(:, 1), win, [], [], Fs);
+  [px_ge, ~] = pwelch(m_e(:, 2), win, [], [], Fs);
+#+end_src
+
+
+The results are shown on figure [[fig:floor_marble_psd_z]] for the Z direction, figure [[fig:floor_marble_psd_n]] for the north direction, and figure [[fig:floor_marble_psd_e]] for the east direction.
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(f, sqrt(px_fz), 'DisplayName', 'Floor - Z');
+  plot(f, sqrt(px_gz), 'DisplayName', 'Granite - Z');
+  hold off;
+  set(gca, 'xscale', 'log');
+  set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  xlim([0.1, 500]);
+#+end_src
+
+#+NAME: fig:floor_marble_psd_z
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/floor_marble_psd_z.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:floor_marble_psd_z
+#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - Z direction
+#+RESULTS: fig:floor_marble_psd_z
+[[file:figs/floor_marble_psd_z.png]]
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(f, sqrt(px_fn), 'DisplayName', 'Floor - N');
+  plot(f, sqrt(px_gn), 'DisplayName', 'Granite - N');
+  hold off;
+  set(gca, 'xscale', 'log');
+  set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  xlim([0.1, 500]);
+#+end_src
+
+#+NAME: fig:floor_marble_psd_n
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/floor_marble_psd_n.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:floor_marble_psd_n
+#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - N direction
+#+RESULTS: fig:floor_marble_psd_n
+[[file:figs/floor_marble_psd_n.png]]
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(f, sqrt(px_fe), 'DisplayName', 'Floor - E');
+  plot(f, sqrt(px_ge), 'DisplayName', 'Granite - E');
+  hold off;
+  set(gca, 'xscale', 'log');
+  set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  xlim([0.1, 500]);
+#+end_src
+
+#+NAME: fig:floor_marble_psd_e
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/floor_marble_psd_e.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:floor_marble_psd_e
+#+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - E direction
+#+RESULTS: fig:floor_marble_psd_e
+[[file:figs/floor_marble_psd_e.png]]
+
+** Compute the transfer function from floor motion to ground motion
+We now compute the transfer function from the floor motion to the granite motion.
+
+The result is shown on figure [[fig:tf_granite]].
+#+begin_src matlab :results none
+  [TZ, f] = tfestimate(m_z(:, 1), -m_z(:, 2), win, [], [], Fs);
+  [TN, ~] = tfestimate(m_n(:, 1), -m_n(:, 2), win, [], [], Fs);
+  [TE, ~] = tfestimate(m_e(:, 1), -m_e(:, 2), win, [], [], Fs);
+#+end_src
+
+#+begin_src matlab :results none :exports none
+  figure;
+  ax1 = subplot(2, 1, 1);
+  hold on;
+  plot(f, abs(TZ), 'DisplayName', 'Z');
+  plot(f, abs(TN), 'DisplayName', 'N');
+  plot(f, abs(TE), 'DisplayName', 'E');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  set(gca, 'XTickLabel',[]);
+  ylabel('Magnitude');
+  legend('Location', 'southwest');
+
+  ax2 = subplot(2, 1, 2);
+  hold on;
+  plot(f, mod(180+180/pi*phase(TZ), 360)-180);
+  plot(f, mod(180+180/pi*phase(TN), 360)-180);
+  plot(f, mod(180+180/pi*phase(TE), 360)-180);
+  hold off;
+  set(gca, 'xscale', 'log');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([10, 100]);
+#+end_src
+
+#+NAME: fig:tf_granite
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/tf_granite.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:tf_granite
+#+CAPTION: Transfer function from floor motion to granite motion
+#+RESULTS: fig:tf_granite
+[[file:figs/tf_granite.png]]
+
+** Conclusion
+#+begin_important
+  - We see resonance of the granite at 33Hz in the horizontal directions
+  - We see two resonances for the z direction: at 60Hz and 75Hz
+#+end_important

dynamical-meas-granite/matlab/marble_dynamics.m

@@ -0,0 +1,153 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+% Load data
+
+m_z = load('mat/data_037.mat', 'data'); m_z = m_z.data;
+m_n = load('mat/data_038.mat', 'data'); m_n = m_n.data;
+m_e = load('mat/data_039.mat', 'data'); m_e = m_e.data;
+
+% Time domain plots
+
+figure;
+subplot(1, 3, 1);
+hold on;
+plot(m_z(:, 3), m_z(:, 2), 'DisplayName', 'Marble - Z');
+plot(m_z(:, 3), m_z(:, 1), 'DisplayName', 'Floor - Z');
+hold off;
+xlabel('Time [s]'); ylabel('Voltage [V]');
+xlim([0, 100]); ylim([-2 2]);
+legend('Location', 'northeast');
+
+subplot(1, 3, 2);
+hold on;
+plot(m_n(:, 3), m_n(:, 2), 'DisplayName', 'Marble - N');
+plot(m_n(:, 3), m_n(:, 1), 'DisplayName', 'Floor - N');
+hold off;
+xlabel('Time [s]'); ylabel('Voltage [V]');
+xlim([0, 100]); ylim([-2 2]);
+legend('Location', 'northeast');
+
+subplot(1, 3, 3);
+hold on;
+plot(m_e(:, 3), m_e(:, 2), 'DisplayName', 'Marble - E');
+plot(m_e(:, 3), m_e(:, 1), 'DisplayName', 'Floor - E');
+hold off;
+xlabel('Time [s]'); ylabel('Voltage [V]');
+xlim([0, 100]); ylim([-2 2]);
+legend('Location', 'northeast');
+
+% Compute the power spectral densities
+% We first compute some parameters that will be used for the PSD computation.
+
+dt = m_z(2, 3)-m_z(1, 3);
+
+Fs = 1/dt; % [Hz]
+
+win = hanning(ceil(10*Fs));
+
+
+
+% Then we compute the Power Spectral Density using =pwelch= function.
+
+[px_fz, f] = pwelch(m_z(:, 1), win, [], [], Fs);
+[px_gz, ~] = pwelch(m_z(:, 2), win, [], [], Fs);
+
+[px_fn, ~] = pwelch(m_n(:, 1), win, [], [], Fs);
+[px_gn, ~] = pwelch(m_n(:, 2), win, [], [], Fs);
+
+[px_fe, ~] = pwelch(m_e(:, 1), win, [], [], Fs);
+[px_ge, ~] = pwelch(m_e(:, 2), win, [], [], Fs);
+
+
+
+
+% The results are shown on figure [[fig:floor_marble_psd_z]] for the Z direction, figure [[fig:floor_marble_psd_n]] for the north direction, and figure [[fig:floor_marble_psd_e]] for the east direction.
+
+
+figure;
+hold on;
+plot(f, sqrt(px_fz), 'DisplayName', 'Floor - Z');
+plot(f, sqrt(px_gz), 'DisplayName', 'Granite - Z');
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+legend('Location', 'southwest');
+xlim([0.1, 500]);
+
+
+
+% #+NAME: fig:floor_marble_psd_z
+% #+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - Z direction
+% #+RESULTS: fig:floor_marble_psd_z
+% [[file:figs/floor_marble_psd_z.png]]
+
+
+figure;
+hold on;
+plot(f, sqrt(px_fn), 'DisplayName', 'Floor - N');
+plot(f, sqrt(px_gn), 'DisplayName', 'Granite - N');
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+legend('Location', 'southwest');
+xlim([0.1, 500]);
+
+
+
+% #+NAME: fig:floor_marble_psd_n
+% #+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - N direction
+% #+RESULTS: fig:floor_marble_psd_n
+% [[file:figs/floor_marble_psd_n.png]]
+
+
+figure;
+hold on;
+plot(f, sqrt(px_fe), 'DisplayName', 'Floor - E');
+plot(f, sqrt(px_ge), 'DisplayName', 'Granite - E');
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
+legend('Location', 'southwest');
+xlim([0.1, 500]);
+
+% Compute the transfer function from floor motion to ground motion
+% We now compute the transfer function from the floor motion to the granite motion.
+
+% The result is shown on figure [[fig:tf_granite]].
+
+[TZ, f] = tfestimate(m_z(:, 1), -m_z(:, 2), win, [], [], Fs);
+[TN, ~] = tfestimate(m_n(:, 1), -m_n(:, 2), win, [], [], Fs);
+[TE, ~] = tfestimate(m_e(:, 1), -m_e(:, 2), win, [], [], Fs);
+
+figure;
+ax1 = subplot(2, 1, 1);
+hold on;
+plot(f, abs(TZ), 'DisplayName', 'Z');
+plot(f, abs(TN), 'DisplayName', 'N');
+plot(f, abs(TE), 'DisplayName', 'E');
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('Magnitude');
+legend('Location', 'southwest');
+
+ax2 = subplot(2, 1, 2);
+hold on;
+plot(f, mod(180+180/pi*phase(TZ), 360)-180);
+plot(f, mod(180+180/pi*phase(TN), 360)-180);
+plot(f, mod(180+180/pi*phase(TE), 360)-180);
+hold off;
+set(gca, 'xscale', 'log');
+ylim([-180, 180]);
+yticks([-180, -90, 0, 90, 180]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+
+linkaxes([ax1,ax2],'x');
+xlim([10, 100]);