#+TITLE: Encoder - Test Bench
:DRAWER:
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+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:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :tangle script.m
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:

* Experimental Setup
The experimental Setup is schematically represented in Figure [[fig:exp_setup_schematic]].

The mass can be vertically moved using the amplified piezoelectric actuator.
The displacement of the mass (relative to the mechanical frame) is measured both by the interferometer and by the encoder.

#+name: fig:exp_setup_schematic
#+caption: Schematic of the Experiment
[[file:figs/exp_setup_schematic.png]]

#+name: fig:encoder_side_view
#+caption: Side View of the encoder
[[file:figs/IMG_20201023_153905.jpg]]

#+name: fig:encoder_front_view
#+caption: Front View of the encoder
[[file:figs/IMG_20201023_153914.jpg]]

* Huddle Test
** Introduction                                                      :ignore:
The goal in this section is the estimate the noise of both the encoder and the intereferometer.

** 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
  load('mat/int_enc_huddle_test.mat', 'interferometer', 'encoder', 't');
#+end_src

#+begin_src matlab
  interferometer = detrend(interferometer, 0);
  encoder = detrend(encoder, 0);
#+end_src

** Time Domain Results
#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, encoder, 'DisplayName', 'Encoder')
  plot(t, interferometer, 'DisplayName', 'Interferometer')
  hold off;
  xlabel('Time [s]'); ylabel('Displacement [m]');
  legend('location', 'northeast');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/huddle_test_time_domain.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:huddle_test_time_domain
#+caption: Huddle test - Time domain signals
#+RESULTS:
[[file:figs/huddle_test_time_domain.png]]

#+begin_src matlab
  G_lpf = 1/(1 + s/2/pi/10);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, lsim(G_lpf, encoder, t), 'DisplayName', 'Encoder')
  plot(t, lsim(G_lpf, interferometer, t), 'DisplayName', 'Interferometer')
  hold off;
  xlabel('Time [s]'); ylabel('Displacement [m]');
  legend('location', 'northeast');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/huddle_test_time_domain_filtered.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:huddle_test_time_domain_filtered
#+caption: Huddle test - Time domain signals filtered with a LPF at 10Hz
#+RESULTS:
[[file:figs/huddle_test_time_domain_filtered.png]]

** Frequency Domain Noise
#+begin_src matlab
  Ts = 1e-4;
  win = hann(ceil(10/Ts));

  [p_i, f] = pwelch(interferometer, win, [], [], 1/Ts);
  [p_e, ~] = pwelch(encoder,        win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(p_i), 'DisplayName', 'Interferometer');
  plot(f, sqrt(p_e), 'DisplayName', 'Encoder');
  hold off;
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
  legend();
  xlim([1e-1, 5e3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/huddle_test_asd.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:huddle_test_asd
#+caption: Amplitude Spectral Density of the signals during the Huddle test
#+RESULTS:
[[file:figs/huddle_test_asd.png]]

* Comparison Interferometer / Encoder
** Introduction                                                      :ignore:
The goal here is to make sure that the interferometer and encoder measurements are coherent.
We may see non-linearity in the interferometric measurement.

** 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
  load('mat/int_enc_comp.mat', 'interferometer', 'encoder', 'u', 't');
#+end_src

#+begin_src matlab
  interferometer = detrend(interferometer, 0);
  encoder = detrend(encoder, 0);
  u = detrend(u, 0);
#+end_src

** Time Domain Results
#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, encoder,        '-',  'DisplayName', 'Encoder')
  plot(t, interferometer, '--', 'DisplayName', 'Interferometer')
  hold off;
  xlabel('Time [s]'); ylabel('Displacement [m]');
  legend('location', 'northeast');
  xlim([50, 52])
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/int_enc_one_cycle.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:int_enc_one_cycle
#+caption: One cycle measurement
#+RESULTS:
[[file:figs/int_enc_one_cycle.png]]

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t, encoder - interferometer, 'DisplayName', 'Difference')
  hold off;
  xlabel('Time [s]'); ylabel('Displacement [m]');
  legend('location', 'northeast');
  xlim([50, 52])
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/int_enc_one_cycle_error.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:int_enc_one_cycle_error
#+caption: Difference between the Encoder and the interferometer during one cycle
#+RESULTS:
[[file:figs/int_enc_one_cycle_error.png]]

** Difference between Encoder and Interferometer as a function of time
#+begin_src matlab
  Ts = 1e-4;
  d_i_mean = reshape(interferometer, [2/Ts floor(Ts/2*length(interferometer))]);
  d_e_mean = reshape(encoder,        [2/Ts floor(Ts/2*length(encoder))]);
#+end_src

#+begin_src matlab
  w0 = 2*pi*5; % [rad/s]
  xi = 0.7;

  G_lpf = 1/(1 + 2*xi/w0*s + s^2/w0^2);

  d_err_mean = reshape(lsim(G_lpf, encoder - interferometer, t), [2/Ts floor(Ts/2*length(encoder))]);
  d_err_mean = d_err_mean - mean(d_err_mean);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  for i_i = 1:size(d_err_mean, 2)
    plot(t(1:size(d_err_mean, 1)), d_err_mean(:, i_i), 'k-')
  end
  plot(t(1:size(d_err_mean, 1)), mean(d_err_mean, 2), 'r-')
  hold off;
  xlabel('Time [s]'); ylabel('Displacement [m]');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/int_enc_error_mean_time.pdf', 'width', 'wide', 'height', 'normal', 'pdf', false);
#+end_src

#+name: fig:int_enc_error_mean_time
#+caption: Difference between the two measurement in the time domain, averaged for all the cycles
#+RESULTS:
[[file:figs/int_enc_error_mean_time.png]]

** Difference between Encoder and Interferometer as a function of position
Compute the mean of the interferometer measurement corresponding to each of the encoder measurement.

#+begin_src matlab
  [e_sorted, ~, e_ind] = unique(encoder);

  i_mean = zeros(length(e_sorted), 1);
  for i = 1:length(e_sorted)
    i_mean(i) = mean(interferometer(e_ind == i));
  end

  i_mean_error = (i_mean - e_sorted);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  % plot(encoder, interferometer - encoder, 'k.', 'DisplayName', 'Difference')
  plot(1e6*(e_sorted), 1e9*(i_mean_error))
  hold off;
  xlabel('Encoder Measurement [$\mu m$]'); ylabel('Measrement Error [nm]');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/int_enc_error_mean_position.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

#+name: fig:int_enc_error_mean_position
#+caption: Difference between the two measurement as a function of the measured position by the encoder, averaged for all the cycles
#+RESULTS:
[[file:figs/int_enc_error_mean_position.png]]

The period of the non-linearity seems to be $1.53 \mu m$ which corresponds to the wavelength of the Laser.

#+begin_src matlab
  win_length = 1530; % length of the windows (corresponds to 1.53 um)
  num_avg = floor(length(e_sorted)/win_length); % number of averaging

  i_init = ceil((length(e_sorted) - win_length*num_avg)/2); % does not start at the extremity

  e_sorted_mean_over_period = mean(reshape(i_mean_error(i_init:i_init+win_length*num_avg-1), [win_length num_avg]), 2);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(1e-3*(0:win_length-1), e_sorted_mean_over_period)
  hold off;
  xlabel('Displacement [$\mu m$]'); ylabel('Measurement Non-Linearity [nm]');
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/int_non_linearity_period_wavelength.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:int_non_linearity_period_wavelength
#+caption: Non-Linearity of the Interferometer over the period of the wavelength
#+RESULTS:
[[file:figs/int_non_linearity_period_wavelength.png]]

* Identification
** 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
  load('mat/int_enc_id_noise_bis.mat', 'interferometer', 'encoder', 'u', 't');
#+end_src

#+begin_src matlab
  interferometer = detrend(interferometer, 0);
  encoder = detrend(encoder, 0);
  u = detrend(u, 0);
#+end_src

** Identification
#+begin_src matlab
  Ts = 1e-4; % Sampling Time [s]
  win = hann(ceil(10/Ts));
#+end_src

#+begin_src matlab
  [tf_i_est, f] = tfestimate(u, interferometer, win, [], [], 1/Ts);
  [co_i_est, ~] = mscohere(u, interferometer, win, [], [], 1/Ts);

  [tf_e_est, ~] = tfestimate(u, encoder, win, [], [], 1/Ts);
  [co_e_est, ~] = mscohere(u, encoder, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, co_i_est, '-')
  plot(f, co_e_est, '-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Coherence'); xlabel('Frequency [Hz]');
  hold off;
  xlim([0.5, 5e3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/identification_dynamics_coherence.pdf', 'width', 'normal', 'height', 'normal');
#+end_src

#+name: fig:identification_dynamics_coherence
#+caption:
#+RESULTS:
[[file:figs/identification_dynamics_coherence.png]]


#+begin_src matlab :exports none
  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  ax1 = nexttile;
  hold on;
  plot(f, abs(tf_i_est), '-', 'DisplayName', 'Int')
  plot(f, abs(tf_e_est), '-', 'DisplayName', 'Enc')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;
  ylim([1e-7, 3e-4]);

  ax2 = nexttile;
  hold on;
  plot(f, 180/pi*angle(tf_i_est), '-')
  plot(f, 180/pi*angle(tf_e_est), '-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;
  yticks(-360:90:360);
  axis padded 'auto x'

  linkaxes([ax1,ax2], 'x');
  xlim([0.5, 5e3]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/identification_dynamics_bode.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:identification_dynamics_bode
#+caption:
#+RESULTS:
[[file:figs/identification_dynamics_bode.png]]