#+TITLE: Test Bench APA95ML
: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:matlab+ :tangle no
#+PROPERTY: header-args:matlab+ :mkdirp yes

#+PROPERTY: header-args:shell  :eval no-export

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

* Introduction                                                        :ignore:

#+name: fig:setup_picture
#+caption: Picture of the Setup
[[file:figs/setup_picture.png]]

#+name: fig:setup_zoom
#+caption: Zoom on the APA
[[file:figs/setup_zoom.png]]

* Setup
:PROPERTIES:
:header-args:matlab+: :tangle matlab/setup_experiment.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** Parameters
#+begin_src matlab
  Ts = 1e-4;
#+end_src

** Filter White Noise
#+begin_src matlab
  Glpf = 1/(1 + s/2/pi/500);

  Gz = c2d(Glpf, Ts, 'tustin');
#+end_src

* Run Experiment and Save Data
:PROPERTIES:
:header-args:matlab+: :tangle matlab/run_experiment.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** Load Data
#+begin_src matlab
  data = SimulinkRealTime.utils.getFileScopeData('data/apa95ml.dat').data;
#+end_src

** Save Data
#+begin_src matlab
  u = data(:, 1); % Input Voltage [V]
  y = data(:, 2); % Output Displacement [m]
  t = data(:, 3); % Time [s]
#+end_src

#+begin_src matlab
  save('./mat/huddle_test.mat', 't', 'u', 'y', 'Glpf');
#+end_src

* Huddle Test
:PROPERTIES:
:header-args:matlab+: :tangle matlab/huddle_test.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** 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                                                       :noexport:
#+begin_src matlab
  load('./mat/huddle_test.mat', 't', 'y');
#+end_src

#+begin_src matlab
  y = y - mean(y(1000:end));
#+end_src

** Time Domain Data
#+begin_src matlab :exports none
  figure;
  plot(t(1000:end), y(1000:end))
  ylabel('Output Displacement [m]'); xlabel('Time [s]');
#+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: Measurement of the Mass displacement during Huddle Test
#+RESULTS:
[[file:figs/huddle_test_time_domain.png]]

** PSD of Measurement Noise
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;
 
  win = hanning(ceil(1*Fs));
#+end_src

#+begin_src matlab
  [pxx, f] = pwelch(y(1000:end), win, [], [], Fs);
#+end_src

#+begin_src matlab :exports none
  figure;
  plot(f, sqrt(pxx));
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
  xlim([1, Fs/2]);
#+end_src

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

#+name: fig:huddle_test_pdf
#+caption: Amplitude Spectral Density of the Displacement during Huddle Test
#+RESULTS:
[[file:figs/huddle_test_pdf.png]]

* Transfer Function Estimation using the DAC as the driver
:PROPERTIES:
:header-args:matlab+: :tangle matlab/tf_estimation.m
:header-args:matlab+: :comments org :mkdirp yes
:END:

** 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                                                       :noexport:
#+begin_src matlab
  ht = load('./mat/huddle_test.mat', 't', 'u', 'y');
  load('./mat/apa95ml_5kg_10V.mat', 't', 'u', 'y');
#+end_src

** Time Domain Data
#+begin_src matlab :exports none
  figure;

  subplot(1,2,1);
  plot(t, u)
  ylabel('Input Voltage [V]'); xlabel('Time [s]');


  subplot(1,2,2);
  plot(t, y)
  ylabel('Output Displacement [m]'); xlabel('Time [s]');
#+end_src

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

#+name: fig:apa95ml_5kg_10V_time_domain
#+caption: Time domain signals during the test
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_time_domain.png]]

** Comparison of the PSD with Huddle Test
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;

  win = hanning(ceil(1*Fs));
#+end_src

#+begin_src matlab
  [pxx, f] = pwelch(y, win, [], [], Fs);
  [pht, ~] = pwelch(ht.y, win, [], [], Fs);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(pxx), 'DisplayName', '5kg');
  plot(f, sqrt(pht), 'DisplayName', 'Huddle Test');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
  legend('location', 'northeast');
  xlim([1, Fs/2]);
#+end_src

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

#+name: fig:apa95ml_5kg_10V_pdf_comp_huddle
#+caption: Comparison of the ASD for the identification test and the huddle test
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_pdf_comp_huddle.png]]

** Compute TF estimate and Coherence
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;
#+end_src

#+begin_src matlab
  win = hann(ceil(1/Ts));

  [tf_est, f] = tfestimate(u, -y, win, [], [], 1/Ts);
  [co_est, ~] = mscohere(  u, -y, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;

  hold on;
  plot(f, co_est, 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Coherence'); xlabel('Frequency [Hz]');
  hold off;
  xlim([10, 5e3]);
#+end_src

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

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

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*angle(tf_est), 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;

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

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

#+name: fig:apa95ml_5kg_10V_tf
#+caption: Estimation of the transfer function from input voltage to displacement
#+RESULTS:
[[file:figs/apa95ml_5kg_10V_tf.png]]

** Comparison with the FEM model
#+begin_src matlab
  load('mat/fem_model_5kg.mat', 'Ghm');
#+end_src

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'DisplayName', 'Identification')
  plot(f, abs(squeeze(freqresp(Ghm, f, 'Hz'))), 'DisplayName', 'FEM')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend('location', 'northeast')
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*angle(tf_est))
  plot(f, 180/pi*angle(squeeze(freqresp(Ghm, f, 'Hz'))))
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;

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

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

#+name: fig:apa95ml_5kg_comp_fem
#+caption: Comparison of the identified transfer function and the one estimated from the FE model
#+RESULTS:
[[file:figs/apa95ml_5kg_comp_fem.png]]

** Conclusion                                                        :ignore:
#+begin_important
The problem comes from the fact that the piezo is driven directly by the DAC that cannot deliver enought current.
In the next section, a current amplifier is used.
#+end_important

* Transfer Function Estimation using the PI Amplifier
** 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                                                       :noexport:
#+begin_src matlab
  ht = load('./mat/huddle_test.mat', 't', 'u', 'y');
  load('./mat/apa95ml_5kg_Amp_E505.mat', 't', 'u', 'um', 'y');
#+end_src

#+begin_src matlab
  u  = u  - mean(u);
  um = um - mean(um);
  y  = y  - mean(y);

  ht.u  = ht.u  - mean(ht.u);
  ht.y  = ht.y  - mean(ht.y);
#+end_src

** Comparison of the PSD with Huddle Test
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;

  win = hanning(ceil(1*Fs));
#+end_src

#+begin_src matlab
  [pxx, f] = pwelch(y, win, [], [], Fs);
  [pht, ~] = pwelch(ht.y, win, [], [], Fs);
#+end_src

#+begin_src matlab :exports none
  figure;
  hold on;
  plot(f, sqrt(pxx), 'DisplayName', '5kg');
  plot(f, sqrt(pht), 'DisplayName', 'Huddle Test');
  hold off;
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
  legend('location', 'southwest');
  xlim([1, Fs/2]);
#+end_src

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

#+name: fig:apa95ml_5kg_PI_pdf_comp_huddle
#+caption: Comparison of the ASD for the identification test and the huddle test
#+RESULTS:
[[file:figs/apa95ml_5kg_PI_pdf_comp_huddle.png]]

** Compute TF estimate and Coherence
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;
#+end_src

#+begin_src matlab
  win = hann(ceil(1/Ts));

  [tf_est, f] = tfestimate(u,  -y, win, [], [], 1/Ts);
  [tf_um , ~] = tfestimate(um, -y, win, [], [], 1/Ts);
  [co_est, ~] = mscohere(  um, -y, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;

  hold on;
  plot(f, co_est, 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Coherence'); xlabel('Frequency [Hz]');
  hold off;
  xlim([10, 5e3]);
#+end_src

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

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

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'DisplayName', 'Input Voltage')
  plot(f, abs(tf_um), 'DisplayName', 'Monitor Voltage')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  hold off;
  legend('location', 'southwest')

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*unwrap(angle(tf_est)))
  plot(f, 180/pi*unwrap(angle(tf_um)))
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;
  ylim([-540, 0]);
  yticks(-540:90:0);

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

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

#+name: fig:apa95ml_5kg_PI_tf
#+caption: Estimation of the transfer function from input voltage to displacement
#+RESULTS:
[[file:figs/apa95ml_5kg_PI_tf.png]]

** Comparison with the FEM model
#+begin_src matlab
  load('mat/fem_model_5kg.mat', 'Ghm');
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(0, 4, 1000);
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, 1/10*170/1400*abs(tf_um), 'DisplayName', 'Identification')
  plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), 'DisplayName', 'FEM')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend('location', 'northeast')
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f,     180/pi*unwrap(angle(tf_um)))
  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))))
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;
  ylim([-540, 0]);
  yticks(-540:90:0);

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

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

#+name: fig:apa95ml_5kg_pi_comp_fem
#+caption: Comparison of the identified transfer function and the one estimated from the FE model
#+RESULTS:
[[file:figs/apa95ml_5kg_pi_comp_fem.png]]
* Transfer function of the PI Amplifier
** 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                                                       :noexport:
#+begin_src matlab
  load('./mat/apa95ml_5kg_Amp_E505.mat', 't', 'u', 'um');
#+end_src

** Compute TF estimate and Coherence
#+begin_src matlab
  Ts = t(end)/(length(t)-1);
  Fs = 1/Ts;
#+end_src

The coherence and the transfer function are estimate from the voltage input of the PI amplifier to its voltage inputs.

The coherence is very good as expected (Figure [[fig:PI_E505_coh]]).

The transfer function show a low pass filter behavior with a lot of phase drop (Figure [[fig:PI_E505_tf]]).

#+begin_src matlab
  win = hann(ceil(10/Ts));

  [tf_est, f] = tfestimate(u, um, win, [], [], 1/Ts);
  [co_est, ~] = mscohere(  u, um, win, [], [], 1/Ts);
#+end_src

#+begin_src matlab :exports none
  figure;

  hold on;
  plot(f, co_est, 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
  ylabel('Coherence'); xlabel('Frequency [Hz]');
  hold off;
  xlim([10, 5e3]);
#+end_src

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

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

#+begin_src matlab :exports none
  figure;
  ax1 = subplot(2, 1, 1);
  hold on;
  plot(f, abs(tf_est), 'k-')
  set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  hold off;

  ax2 = subplot(2, 1, 2);
  hold on;
  plot(f, 180/pi*angle(tf_est), 'k-')
  set(gca, 'Xscale', 'lin'); set(gca, 'Yscale', 'lin');
  ylabel('Phase'); xlabel('Frequency [Hz]');
  hold off;
  ylim([-180, 180]);
  yticks(-180:90:180);

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

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

#+name: fig:PI_E505_tf
#+caption: Estimation of the transfer function from input voltage to displacement
#+RESULTS:
[[file:figs/PI_E505_tf.png]]

The delay can be estimated as follow (in ms):
#+begin_src matlab :results replace value
  finddelay(u, um)*(1000*Ts)
#+end_src

#+RESULTS:
: 0.4

This most probably corresponds to a FIR filter.