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

Nothing is then to the actuator such that the relative motion between the mass and the frame is as small as possible.
Ideally, a mechanical part would clamp the two together, we here suppose that the APA is still enough to clamp the two together.

** 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), 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]]

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

* Change of Stiffness due to Sensors stack being open/closed circuit
** 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
  oc = load('./mat/identification_open_circuit.mat', 't', 'encoder', 'u');
  sc = load('./mat/identification_short_circuit.mat', 't', 'encoder', 'u');
#+end_src

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

#+begin_src matlab
  [tf_oc_est, f] = tfestimate(oc.u, oc.encoder, win, [], [], 1/Ts);
  [co_oc_est, ~] = mscohere(  oc.u, oc.encoder, win, [], [], 1/Ts);

  [tf_sc_est, ~] = tfestimate(sc.u, sc.encoder, win, [], [], 1/Ts);
  [co_sc_est, ~] = mscohere(  sc.u, sc.encoder, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, co_oc_est, '-')
  plot(f, co_sc_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/stiffness_force_sensor_coherence.pdf', 'width', 'wide', 'height', 'normal');
#+end_src

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


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

  ax1 = nexttile;
  hold on;
  plot(f, abs(tf_oc_est), '-', 'DisplayName', 'Open-Circuit')
  plot(f, abs(tf_sc_est), '-', 'DisplayName', 'Short-Circuit')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
  hold off;
  ylim([1e-7, 3e-4]);
  legend('location', 'southwest');

  ax2 = nexttile;
  hold on;
  plot(f, 180/pi*angle(tf_oc_est), '-')
  plot(f, 180/pi*angle(tf_sc_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/stiffness_force_sensor_bode.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

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

#+begin_src matlab :tangle no :exports results :results file replace
  xlim([180, 280]);
  exportFig('figs/stiffness_force_sensor_bode_zoom.pdf', 'width', 'small', 'height', 'tall');
#+end_src

#+name: fig:stiffness_force_sensor_bode_zoom
#+caption: Zoom on the change of resonance
#+RESULTS:
[[file:figs/stiffness_force_sensor_bode_zoom.png]]

#+begin_important
  The change of resonance frequency / stiffness is very small and is not important here.
#+end_important

* Generated Number of Charge / Voltage
** Introduction                                                      :ignore:
Two stacks are used as actuator (in parallel) and one stack is used as sensor.

The amplifier gain is 20V/V (Cedrat LA75B).

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

** Steps
#+begin_src matlab
  load('./mat/force_sensor_steps.mat', 't', 'encoder', 'u', 'v');
#+end_src

#+begin_src matlab
  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
  nexttile;
  plot(t, v);
  xlabel('Time [s]'); ylabel('Measured voltage [V]');
  nexttile;
  plot(t, u);
  xlabel('Time [s]'); ylabel('Actuator Voltage [V]');
#+end_src

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

#+name: fig:force_sen_steps_time_domain
#+caption: Time domain signal during the 3 actuator voltage steps
#+RESULTS:
[[file:figs/force_sen_steps_time_domain.png]]

Three steps are performed at the following time intervals:
- 2.5, 23
- 23.8, 35
- 35.8, 50

Fit function:
#+begin_src matlab
  f = @(b,x) b(1).*exp(b(2).*x) + b(3);
#+end_src

We are interested by the =b(2)= term, which is the time constant of the exponential.

#+begin_src matlab
  tau = zeros(3,1);
  V0  = zeros(3,1);
#+end_src

#+begin_src matlab
  t_cur = t(2.5 < t & t < 23);
  t_cur = t_cur - t_cur(1);
  y_cur = v(2.5 < t & t < 23);

  nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
  B0 = [0.5, -0.15, 2.2];                        % Choose Appropriate Initial Estimates
  [B,rnrm] = fminsearch(nrmrsd, B0);     % Estimate Parameters ‘B’

  tau(1) = 1/B(2);
  V0(1)  = B(3);
#+end_src

#+begin_src matlab
  t_cur = t(23.8 < t & t < 35);
  t_cur = t_cur - t_cur(1);
  y_cur = v(23.8 < t & t < 35);

  nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
  B0 = [0.5, -0.15, 2.2];                 % Choose Appropriate Initial Estimates
  [B,rnrm] = fminsearch(nrmrsd, B0);      % Estimate Parameters ‘B’

  tau(2) = 1/B(2);
  V0(2)  = B(3);
#+end_src

#+begin_src matlab
  t_cur = t(35.8 < t);
  t_cur = t_cur - t_cur(1);
  y_cur = v(35.8 < t);

  nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
  B0 = [0.5, -0.15, 2.2];                 % Choose Appropriate Initial Estimates
  [B,rnrm] = fminsearch(nrmrsd, B0);     % Estimate Parameters ‘B’

  tau(3) = 1/B(2);
  V0(3)  = B(3);
#+end_src

#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
  data2orgtable([abs(tau), V0], {}, {'$tau$ [s]', '$V_0$ [V]'}, ' %.2f ');
#+end_src

#+RESULTS:
| $tau$ [s] | $V_0$ [V] |
|-----------+-----------|
|      6.47 |      2.26 |
|      6.76 |      2.26 |
|      6.49 |      2.25 |

With the capacitance being $C = 4.4 \mu F$, the internal impedance of the Speedgoat ADC can be computed as follows:
#+begin_src matlab
  Cp = 4.4e-6; % [F]
  Rin = abs(mean(tau))/Cp;
#+end_src

#+begin_src matlab :results value replace :exports results
  ans = Rin
#+end_src

#+RESULTS:
: 1494100.0

The input impedance of the Speedgoat's ADC should then be close to $1.5\,M\Omega$ (specified at $1\,M\Omega$).

#+begin_important
  How can we explain the voltage offset?
#+end_important

As shown in Figure [[fig:piezo_sensor_model_instrumentation]] (taken from cite:reza06_piezoel_trans_vibrat_contr_dampin), an input voltage offset is due to the input bias current $i_n$.

#+name: fig:piezo_sensor_model_instrumentation
#+caption: Model of a piezoelectric transducer (left) and instrumentation amplifier (right)
[[file:figs/piezo_sensor_model_instrumentation.png]]

The estimated input bias current is then:
#+begin_src matlab
  in = mean(V0)/Rin;
#+end_src

#+begin_src matlab :results value replace :exports results
   ans = in
#+end_src

#+RESULTS:
: 1.5119e-06

An additional resistor in parallel with $R_{in}$ would have two effects:
- reduce the input voltage offset
  \[ V_{off} = \frac{R_a R_{in}}{R_a + R_{in}} i_n \]
- increase the high pass corner frequency $f_c$
  \[ C_p \frac{R_{in}R_a}{R_{in} + R_a} = \tau_c = \frac{1}{f_c} \]
  \[ R_a = \frac{R_i}{f_c C_p R_i - 1} \]


If we allow the high pass corner frequency to be equals to 3Hz:
#+begin_src matlab
  fc = 3;
  Ra = Rin/(fc*C*Rin - 1);
#+end_src

#+begin_src matlab :results value replace :exports results
  ans = Ra
#+end_src

#+RESULTS:
: 79804

With this parallel resistance value, the voltage offset would be:
#+begin_src matlab
  V_offset = Ra*Rin/(Ra + Rin) * in;
#+end_src

#+begin_src matlab :results value replace :exports results
  ans = V_offset
#+end_src

#+RESULTS:
: 0.11454

Which is much more acceptable.

** Sinus
#+begin_src matlab
  load('./mat/force_sensor_sin.mat', 't', 'encoder', 'u', 'v');

  u       = u(t>25);
  v       = v(t>25);
  encoder = encoder(t>25) - mean(encoder(t>25));
  t       = t(t>25);
#+end_src

The driving voltage is a sinus at 0.5Hz centered on 3V and with an amplitude of 3V (Figure [[fig:force_sensor_sin_u]]).

#+begin_src matlab :exports none
  figure;
  plot(t, u)
  xlabel('Time [s]'); ylabel('Control Voltage [V]');
#+end_src

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

#+name: fig:force_sensor_sin_u
#+caption: Driving Voltage
#+RESULTS:
[[file:figs/force_sensor_sin_u.png]]

The full stroke as measured by the encoder is:
#+begin_src matlab :results value replace
  max(encoder)-min(encoder)
#+end_src

#+RESULTS:
: 5.005e-05

Its signal is shown in Figure [[fig:force_sensor_sin_encoder]].

#+begin_src matlab :exports none
  figure;
  plot(t, encoder)
  xlabel('Time [s]'); ylabel('Encoder [m]');
#+end_src

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

#+name: fig:force_sensor_sin_encoder
#+caption: Encoder measurement
#+RESULTS:
[[file:figs/force_sensor_sin_encoder.png]]

The generated voltage by the stack is shown in Figure

#+begin_src matlab :exports none
  figure;
  plot(t, v)
  xlabel('Time [s]'); ylabel('Force Sensor Output [V]');
#+end_src

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

#+name: fig:force_sensor_sin_stack
#+caption: Voltage measured on the stack used as a sensor
#+RESULTS:
[[file:figs/force_sensor_sin_stack.png]]

The capacitance of the stack is
#+begin_src matlab
  Cp = 4.4e-6; % [F]
#+end_src

The corresponding generated charge is then shown in Figure [[fig:force_sensor_sin_charge]].
#+begin_src matlab :exports none
  figure;
  plot(t, 1e6*Cp*(v-mean(v)))
  xlabel('Time [s]'); ylabel('Generated Charge [$\mu C$]');
#+end_src

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

#+name: fig:force_sensor_sin_charge
#+caption: Generated Charge
#+RESULTS:
[[file:figs/force_sensor_sin_charge.png]]


The relation between the generated voltage and the measured displacement is almost linear as shown in Figure [[fig:force_sensor_linear_relation]].

#+begin_src matlab
  b1 = encoder\(v-mean(v));
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(encoder, v-mean(v), 'DisplayName', 'Measured Voltage');
  plot(encoder, encoder*b1, 'DisplayName', sprintf('Linear Fit: $U_s \\approx %.3f [V/\\mu m] \\cdot d$', 1e-6*abs(b1)));
  hold off;
  xlabel('Measured Displacement [m]'); ylabel('Generated Voltage [V]');
  legend();
#+end_src

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

#+name: fig:force_sensor_linear_relation
#+caption: Almost linear relation between the relative displacement and the generated voltage
#+RESULTS:
[[file:figs/force_sensor_linear_relation.png]]

With a 16bits ADC, the resolution will then be equals to (in [nm]):
#+begin_src matlab :results value replace
  abs((20/2^16)/(b1/1e9))
#+end_src

#+RESULTS:
: 3.9838