#+TITLE: Attocube - Test Bench
:DRAWER:
#+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/custom.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/readtheorg.js"></script>

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :tangle no
#+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:

* Estimation of the Spectral Density of the Attocube Noise
** Introduction                                                      :ignore:

#+name: fig:test-bench-schematic
#+caption: Test Bench Schematic
[[file:figs/test-bench-schematic.png]]

#+name: fig:test-bench-picture
#+caption: Picture of the test bench. The Attocube and mirror are covered by a "bubble sheet"
#+ATTR_ORG: :width 600px
[[file:figs/IMG-7865.JPG]]

** 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

** Long and Slow measurement
The first measurement was made during ~17 hours with a sampling time of $T_s = 0.1\,s$.

#+begin_src matlab
  load('./mat/long_test_plastic.mat', 'x', 't')
  Ts = 0.1; % [s]
#+end_src

#+begin_src matlab :exports none
  figure;
  plot(t/60/60, 1e9*x)
  xlim([0, 17.5]);
  xlabel('Time [h]'); ylabel('Displacement [nm]');
#+end_src

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

#+name: fig:long_meas_time_domain_full
#+caption: Long measurement time domain data
#+RESULTS:
[[file:figs/long_meas_time_domain_full.png]]

Let's fit the data with a step response to a first order low pass filter (Figure [[fig:long_meas_time_domain_fit]]).

#+begin_src matlab
  f = @(b,x) b(1)*(1 - exp(-x/b(2)));

  y_cur = x(t < 17.5*60*60);
  t_cur = t(t < 17.5*60*60);

  nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
  B0 = [400e-9, 2*60*60];                        % Choose Appropriate Initial Estimates
  [B,rnrm] = fminsearch(nrmrsd, B0);     % Estimate Parameters ‘B’
#+end_src

The corresponding time constant is (in [h]):
#+begin_src matlab :results value replace :exports results
  B(2)/60/60
#+end_src

#+RESULTS:
: 2.0658

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(t_cur/60/60, 1e9*y_cur);
  plot(t_cur/60/60, 1e9*f(B, t_cur));
  hold off;
  xlim([0, 17.5])
  xlabel('Time [h]'); ylabel('Displacement [nm]');
#+end_src

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

#+name: fig:long_meas_time_domain_fit
#+caption: Fit of the measurement data with a step response of a first order low pass filter
#+RESULTS:
[[file:figs/long_meas_time_domain_fit.png]]

We can see in Figure [[fig:long_meas_time_domain_full]] that there is a transient period where the measured displacement experiences some drifts.
This is probably due to thermal effects.
We only select the data between =t1= and =t2=.
The obtained displacement is shown in Figure [[fig:long_meas_time_domain_zoom]].

#+begin_src matlab
  t1 = 10.5; t2 = 17.5; % [h]

  x = x(t > t1*60*60 & t < t2*60*60);
  x = x - mean(x);
  t = t(t > t1*60*60 & t < t2*60*60);
  t = t - t(1);
#+end_src

#+begin_src matlab :exports none
  figure;
  plot(t/60/60, 1e9*x);
  xlabel('Time [h]'); ylabel('Measured Displacement [nm]')
#+end_src

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

#+name: fig:long_meas_time_domain_zoom
#+caption: Kept data (removed slow drifts during the first hours)
#+RESULTS:
[[file:figs/long_meas_time_domain_zoom.png]]

The Power Spectral Density of the measured displacement is computed
#+begin_src matlab
  win = hann(ceil(length(x)/20));
  [p_1, f_1] = pwelch(x, win, [], [], 1/Ts);
#+end_src

As a low pass filter was used in the measurement process, we multiply the PSD by the square of the inverse of the filter's norm.
#+begin_src matlab
  G_lpf = 1/(1 + s/2/pi);
  p_1 = p_1./abs(squeeze(freqresp(G_lpf, f_1, 'Hz'))).^2;
#+end_src

Only frequencies below 2Hz are taken into account (high frequency noise will be measured afterwards).
#+begin_src matlab
  p_1 = p_1(f_1 < 2);
  f_1 = f_1(f_1 < 2);
#+end_src

** Short and Fast measurement
An second measurement is done in order to estimate the high frequency noise of the interferometer.
The measurement is done with a sampling time of $T_s = 0.1\,ms$ and a duration of ~100s.

#+begin_src matlab
  load('./mat/short_test_plastic.mat')
  Ts = 1e-4; % [s]
#+end_src

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

The time domain measurement is shown in Figure [[fig:short_meas_time_domain]].

#+begin_src matlab :exports none
  figure;
  plot(t, 1e9*x)
  xlabel('Time [s]'); ylabel('Displacement [nm]');
#+end_src

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

#+name: fig:short_meas_time_domain
#+caption: Time domain measurement with the high sampling rate
#+RESULTS:
[[file:figs/short_meas_time_domain.png]]

The Power Spectral Density of the measured displacement is computed
#+begin_src matlab
  win = hann(ceil(length(x)/20));
  [p_2, f_2] = pwelch(x, win, [], [], 1/Ts);
#+end_src

** Obtained Amplitude Spectral Density of the measured displacement

The computed ASD of the two measurements are combined in Figure [[fig:psd_combined]].

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f_1(8:end), sqrt(p_1(8:end)), 'k-');
  plot(f_2(8:end), sqrt(p_2(8:end)), 'k-');
  hold off;
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
#+end_src

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

#+name: fig:psd_combined
#+caption: Obtained Amplitude Spectral Density of the measured displacement
#+RESULTS:
[[file:figs/psd_combined.png]]

* Effect of the "bubble sheet" and "Aluminium tube"
** Introduction                                                      :ignore:

 #+name: fig:picture-test-bench-aluminium-tube
 #+caption: Aluminium tube used to protect the beam path from disturbances
#+ATTR_ORG: :width 600px
 [[file:figs/IMG-7864.JPG]]

** 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

** Aluminium Tube and Bubble Sheet
#+begin_src matlab
  load('./mat/short_test_plastic.mat');
  Ts = 1e-4; % [s]
#+end_src

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

#+begin_src matlab :exports none
  figure;
  plot(t, 1e9*x)
  xlabel('Time [s]'); ylabel('Displacement [nm]');
#+end_src

#+begin_src matlab
  win = hann(ceil(length(x)/10));
  [p_1, f_1] = pwelch(x, win, [], [], 1/Ts);
#+end_src

** Only Aluminium Tube
#+begin_src matlab
  load('./mat/short_test_alu_tube.mat');
  Ts = 1e-4; % [s]
#+end_src

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

The time domain measurement is shown in Figure [[fig:short_meas_time_domain]].
#+begin_src matlab :exports none
  figure;
  plot(t, 1e9*x)
  xlabel('Time [s]'); ylabel('Displacement [nm]');
#+end_src

#+begin_src matlab
  win = hann(ceil(length(x)/10));
  [p_2, f_2] = pwelch(x, win, [], [], 1/Ts);
#+end_src

** Nothing
#+begin_src matlab
  load('./mat/short_test_without_material.mat');
  Ts = 1e-4; % [s]
#+end_src

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

The time domain measurement is shown in Figure [[fig:short_meas_time_domain]].
#+begin_src matlab :exports none
  figure;
  plot(t, 1e9*x)
  xlabel('Time [s]'); ylabel('Displacement [nm]');
#+end_src

#+begin_src matlab
  win = hann(ceil(length(x)/10));
  [p_3, f_3] = pwelch(x, win, [], [], 1/Ts);
#+end_src

** Comparison
#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f_1(8:end), sqrt(p_1(8:end)), '-', ...
       'DisplayName', 'Alunimium + Bubble');
  plot(f_2(8:end), sqrt(p_2(8:end)), '-', ...
       'DisplayName', 'Aluminium');
  plot(f_3(8:end), sqrt(p_3(8:end)), '-', ...
       'DisplayName', 'nothing');
  hold off;
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
  xlim([1e-1, 5e3]);
  legend('location', 'northeast');
#+end_src

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

#+name: fig:asd_noise_comp_bubble_aluminium
#+caption: Comparison of the noise ASD with and without bubble sheet
#+RESULTS:
[[file:figs/asd_noise_comp_bubble_aluminium.png]]

* Measurement of the Attocube's non-linearity
** Introduction                                                      :ignore:
The measurement setup is shown in Figure [[fig:exp_setup_schematic]].

#+begin_note
Here are the equipment used in the test bench:
- Renishaw Resolution Encoder with 1nm resolution ([[file:doc/L-9517-9448-05-B_Data_sheet_RESOLUTE_BiSS_en.pdf][doc]])
- Attocube interferometer ([[file:doc/IDS3010.pdf][doc]])
- Cedrat Amplified Piezoelectric Actuator APA95ML ([[file:doc/APA95ML.pdf][doc]])
- Voltage Amplifier LA75B ([[file:doc/LA75B.pdf][doc]])
- Speedgoat IO131 with 16bits ADC and DAC ([[file:doc/IO130 IO131 OEM Datasheet.pdf][doc]])
#+end_note

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

A DAC and voltage amplified are used to move the mass with the Amplified Piezoelectric Actuator (APA95ML).
The encoder and the attocube are measure ring the same motion.

As will be shown shortly, this measurement permitted to measure the period non-linearity of the Attocube.

** 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
The measurement data are loaded and the offset are removed using the =detrend= command.

#+begin_src matlab
  load('mat/int_enc_comp.mat', 'interferometer', 'encoder', 'u', 't');
  Ts = 1e-4; % Sampling Time [s]
#+end_src

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

** Time Domain Results
One period of the displacement of the mass as measured by the encoder and interferometer are shown in Figure [[fig:int_enc_one_cycle]].
It consist of the sinusoidal motion at 0.5Hz with an amplitude of approximately $70\mu m$.

The frequency of the motion is chosen such that no resonance in the system is excited.
This should improve the coherence between the measurements made by the encoder and interferometer.

#+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', 'southeast');
  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]]

The difference between the two measurements during the same period is shown in Figure [[fig:int_enc_one_cycle_error]].

#+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
The data is filtered using a second order low pass filter with a cut-off frequency $\omega_0$ as defined below.

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

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

After filtering, the data is "re-shaped" such that we can superimpose all the measured periods as shown in Figure [[fig:int_enc_error_mean_time]].
This gives an idea of the measurement error as given by the Attocube during a $70 \mu m$ motion.
#+begin_src matlab
  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
Figure [[fig:int_enc_error_mean_time]] gives the measurement error as a function of time.
We here wish the compute this measurement error as a function of the position (as measured by the encoer).

To do so, all the attocube measurements corresponding to each position measured by the Encoder (resolution of $1nm$) are averaged.
Figure [[fig:int_enc_error_mean_position]] is obtained where we clearly see an error with a period comparable to the motion range and a much smaller period corresponding to the non-linear period errors that we wish the estimate.
#+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 equal to $765 nm$ which corresponds to half the wavelength of the Laser ($1.53 \mu m$).
For the motion range done here, the non-linearity is measured over ~18 periods which permits to do some averaging.

#+begin_src matlab
  win_length = 1530/2; % length of the windows (corresponds to 765 nm)
  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

The obtained periodic non-linearity is shown in Figure [[fig:int_non_linearity_period_wavelength]].

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(1e-3*(0:win_length-1), 1e9*(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]]