#+TITLE:Measurement of the sample vibrations when rotating the Spindle
:DRAWER:
#+STARTUP: overview

#+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>

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

* Signal Processing
** Matlab Init                                              :noexport:ignore:
#+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

** Pre-processing
#+begin_src matlab :results none
  imax = min([length(data1.t), length(data2.t)]);

  data1.t  = data1.t(1:imax);
  data1.x1 = data1.x1(1:imax);
  data1.x2 = data1.x2(1:imax);

  data2.t  = data2.t(1:imax);
  data2.x1 = data2.x1(1:imax);
  data2.x2 = data2.x2(1:imax);
#+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/dt));
#+end_src

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

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

#+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