#+TITLE: Ground Motion Measurements
:DRAWER:
#+STARTUP: overview

#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html

#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>

#+HTML_MATHJAX: align: center tagside: right font: TeX

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs

#+PROPERTY: header-args:shell  :eval no-export
:END:

* Experimental Setup
The goal here is to compare the ground motion at the location of the micro-station (ID31 beamline at ESRF) with other measurements of the ground motion.

This will also permit to confirm that the device use for the measurement (geophones, amplifiers, ADC) are working well and that the data processing is correct.

One L22 geophone is put on the ID31 floor.
The signal is then filtered with a first order low pass filter with a cut-off frequency of $1kHz$.
Then the signal is amplified by a Voltage Amplifier with the following settings:
- AC/DC option set to DC
- Amplification of 60dB (1000)
- Low pass filter at the output with a cut-off frequency of $1kHz$

On figure [[fig:ground_motion_measurements]] is an overview of multiple measurements made at famous location.

#+name: fig:ground_motion_measurements
#+caption: Power spectral density of ground vibration in the vertical direction at multiple locations ([[https://vibration.desy.de/][source]])
#+attr_html: :width 800px
[[file:./img/ground_motion_measurements.png]]

* Measurement Analysis
  :PROPERTIES:
  :header-args:matlab+: :tangle matlab/ground_meas_id31.m
  :header-args:matlab+: :comments org :mkdirp yes
  :END:
  <<sec:ground_meas_id31>>

** ZIP file containing the data and matlab files                     :ignore:
#+begin_src bash :exports none :results none
  if [ matlab/ground_meas_id31.m -nt data/ground_meas_id31.zip ]; then
    cp matlab/ground_meas_id31.m ground_meas_id31.m;
    zip data/ground_meas_id31 \
        mat/data_028.mat \
        mat/id09_floor_september2018.mat \
        mat/ground_motion_dist.mat \
        ground_meas_id31.m
    rm ground_meas_id31.m;
  fi
#+end_src

#+begin_note
  All the files (data and Matlab scripts) are accessible [[file:data/ground_meas_id31.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
  data = load('mat/data_028.mat', 'data'); data = data.data;
#+end_src

** Time domain plots
#+begin_src matlab
  figure;
  hold on;
  plot(data(:, 3), data(:, 1));
  hold off;
  xlabel('Time [s]'); ylabel('Voltage [V]');
  xlim([0, 100]);
#+end_src

#+NAME: fig:ground_motion_id31_time
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_time.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_id31_time
#+CAPTION: Measurement of the ground motion - Time domain
#+RESULTS: fig:ground_motion_id31_time
[[file:figs/ground_motion_id31_time.png]]

** Computation of the ASD of the measured voltage
#+begin_src matlab :results none
  dt = data(2, 3) - data(1, 3);

  Fs = 1/dt;
  win = hanning(ceil(10*Fs));
#+end_src

#+begin_src matlab :results none
  [px_dc, f] = pwelch(data(:, 1), win, [], [], Fs);
#+end_src

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(px_dc));
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density $\left[\frac{V}{\sqrt{Hz}}\right]$')
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:ground_motion_id31_asd_volt
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_volt.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_id31_asd_volt
#+CAPTION: Amplitude Spectral Density of the measured Voltage
#+RESULTS: fig:ground_motion_id31_asd_volt
[[file:figs/ground_motion_id31_asd_volt.png]]

** Scaling to take into account the sensibility of the geophone and the voltage amplifier
The Geophone used are L22. Their sensibility is shown on figure [[fig:geophone_sensibility]].

#+begin_src matlab :results none
  S0 = 88; % Sensitivity [V/(m/s)]
  f0 = 2; % Cut-off frequency [Hz]

  S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
#+end_src

#+begin_src matlab :results none :exports none
  figure;
  bodeFig({S}, logspace(-1, 2, 1000));
  ylabel('Amplitude $\left[\frac{V}{m/s}\right]$')
#+end_src

#+NAME: fig:geophone_sensibility
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/geophone_sensibility.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:geophone_sensibility
#+CAPTION: Sensibility of the Geophone
#+RESULTS: fig:geophone_sensibility
[[file:figs/geophone_sensibility.png]]

We also take into account the gain of the electronics which is here set to be $60dB$.

#+begin_src matlab :results none
  G0_db = 60; % [dB]

  G0 = 10^(G0_db/20); % [abs]
#+end_src

We divide the ASD measured (in $\text{V}/\sqrt{\text{Hz}}$) by the gain of the voltage amplifier to obtain the ASD of the voltage across the geophone.
We further divide the result by the sensibility of the Geophone to obtain the ASD of the velocity in $m/s/\sqrt{Hz}$.

#+begin_src matlab :results none
  scaling = 1./squeeze(abs(freqresp(G0*S, f, 'Hz')));
#+end_src

** Computation of the ASD of the velocity and displacement
The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd_velocity]].

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(px_dc).*scaling);
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD of the measured Velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:ground_motion_id31_asd_velocity
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_id31_asd_velocity
#+CAPTION: Amplitude Spectral Density of the Velocity
#+RESULTS: fig:ground_motion_id31_asd_velocity
[[file:figs/ground_motion_id31_asd_velocity.png]]

We also plot the ASD in displacement (figure [[fig:ground_motion_id31_asd_displacement]]);

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, (sqrt(px_dc).*scaling)./(2*pi*f));
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:ground_motion_id31_asd_displacement
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_asd_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_id31_asd_displacement
#+CAPTION: Amplitude Spectral Density of the Displacement
#+RESULTS: fig:ground_motion_id31_asd_displacement
[[file:figs/ground_motion_id31_asd_displacement.png]]

And we also plot the PSD of the displacement in $\frac{{\mu u}^2}{Hz}$ as it is a usual unit used (figure [[fig:ground_motion_id31_psd_displacement]]).
One can then compare this curve with the figure [[fig:ground_motion_measurements]].

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f).*1e6).^2);
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('PSD of the measured displacement $\left[\frac{{ \mu m }^2}{Hz}\right]$')
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:ground_motion_id31_psd_displacement
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_id31_psd_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_id31_psd_displacement
#+CAPTION: Power Spectral Density of the measured displacement
#+RESULTS: fig:ground_motion_id31_psd_displacement
[[file:figs/ground_motion_id31_psd_displacement.png]]

** Comparison with other measurements of ground motion
Now we will compare with other measurements.

*** Load the measurement data
First we load the measurement data.
Here we have one measurement of the floor motion made at the ESRF in 2018, and one measurement made at CERN.
#+begin_src matlab
  id09 = load('./mat/id09_floor_september2018.mat');
  cern = load('./mat/ground_motion_dist.mat');
#+end_src

*** Compute PSD of the measurements
We compute the Power Spectral Densities of the measurements.
#+begin_src matlab
  Fs_id09 = 1/(id09.time3(2)-id09.time3(1));
  win_id09 = hanning(ceil(10*Fs_id09));
  [id09_pxx, id09_f] = pwelch(1e-6*id09.x_y_z(:, 3), win_id09, [], [], Fs_id09);
#+end_src

#+begin_src matlab
  Fs_cern = 1/(cern.gm.time(2)-cern.gm.time(1));
  win_cern = hanning(ceil(10*Fs_cern));
  [cern_pxx, cern_f] = pwelch(cern.gm.signal, win_cern, [], [], Fs_cern);
#+end_src

*** Compare PSD of Cern, ID09 and ID31
And we compare all the measurements (figure [[fig:ground_motion_compare]]).

#+begin_src matlab :results silent
  figure;
  hold on;
  plot(id09_f, id09_pxx, 'DisplayName', 'ID09');
  plot(cern_f, cern_pxx, 'DisplayName', 'CERN');
  plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f)).^2, 'k', 'DisplayName', 'ID31');
  plot(f, px_disp);
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('PSD [$m^2/Hz$]');
  legend('Location', 'northeast');
  xlim([0.1, 500]);
#+end_src

#+NAME: fig:ground_motion_compare
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/ground_motion_compare.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:ground_motion_compare
#+CAPTION: Comparison of the PSD of the ground motion measured at different location
#+RESULTS: fig:ground_motion_compare
[[file:figs/ground_motion_compare.png]]