#+TITLE:Measurement of the sample vibrations when rotating the Spindle
: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

* Signal Processing
** 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
Measurement =data_001.mat= corresponds to a measurement where the spindle is not turning and =data_002.mat= where the spindle is turning at 1rpm.
=x1= is the signal coming from the geophone located on the marble and =x2= is the signal from the geophone located on the sample station.

#+begin_src matlab :results none
  data1 = load('mat/data_001.mat', 't', 'x1', 'x2');
  data2 = load('mat/data_002.mat', 't', 'x1', 'x2');
#+end_src

** Time domain Data
#+begin_src matlab :results none
  figure;
  hold on;
  plot(data1.t, data1.x1);
  plot(data2.t, data2.x1);
  hold off;
#+end_src

#+begin_src matlab :results none
  figure;
  hold on;
  plot(data1.t, data1.x2);
  plot(data2.t, data2.x2)
  hold off;
#+end_src

** ASD and Frequency domain data
#+begin_src matlab :results none
  dt = data1.t(2) - data1.t(1);
  Fs = 1/dt;

  windows_psd = hanning(ceil(10*Fs));
#+end_src

#+begin_src matlab :results none
  [pxx1m, f] = pwelch(data1.x1, windows_psd, [], [], Fs); f(1) = []; pxx1m(1) = [];
  [pxx1h, ~] = pwelch(data1.x2, windows_psd, [], [], Fs); pxx1h(1) = [];

  [pxx2m, ~] = pwelch(data2.x1, windows_psd, [], [], Fs); pxx2m(1) = [];
  [pxx2h, ~] = pwelch(data2.x2, windows_psd, [], [], Fs); pxx2h(1) = [];
#+end_src

** Some plots
#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(pxx1m));
  plot(f, sqrt(pxx2m));
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(pxx1h));
  plot(f, sqrt(pxx2h));
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src


#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(pxx2m));
  plot(f, sqrt(pxx2h));
  hold off;
  set(gca, 'xscale', 'log');
  set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, cumtrapz(f, pxx1m))
  plot(f, cumtrapz(f, pxx2m))
  set(gca, 'XScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS [m]')
#+end_src

** 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 frequnecy [Hz]

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

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'))); scaling(1) = 0;
#+end_src

** Computation of the ASD of the velocity
#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, sqrt(pxx1h).*scaling);
  plot(f, sqrt(pxx2h).*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

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, (sqrt(pxx1).*scaling)./(2*pi*f));
  plot(f, (sqrt(pxx2).*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

** RMS value of the difference between the two geophones
We also compute the Power Spectral Density of the difference between the two geophones. This is done in order to estimate the relative displacement of the sample with respect to the granite.
#+begin_src matlab :results none
  [pxxd1, ~] = pwelch(data1.x2-data1.x1, windows_psd, [], [], Fs); pxxd1(1) = [];
  [pxxd2, ~] = pwelch(data2.x2-data2.x1, windows_psd, [], [], Fs); pxxd2(1) = [];
#+end_src

#+begin_src matlab :results none
  figure;
  hold on;
  plot(f, (sqrt(pxxd1).*scaling)./(2*pi*f));
  plot(f, (sqrt(pxxd2).*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

#+begin_src matlab :results none
  psd_d1 = ((sqrt(pxxd1).*scaling)./(2*pi*f)).^2;
  psd_d2 = ((sqrt(pxxd2).*scaling)./(2*pi*f)).^2;

  df = f(2) - f(1);

  figure;
  hold on;
  plot(f, sqrt(cumsum(df.*psd_d1, 'reverse')));
  plot(f, sqrt(cumsum(df.*psd_d2, 'reverse')));
  hold off;
  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
  xlabel('Frequency [Hz]'); ylabel('CAS $\left[m\right]$')
  xlim([0.1, 500]);
#+end_src