From 73dc280ab3406e1e9108356f7f32f87722d4afc9 Mon Sep 17 00:00:00 2001
From: Thomas Dehaeze <dehaeze.thomas@gmail.com>
Date: Thu, 11 Feb 2021 10:09:23 +0100
Subject: [PATCH] Prepare files for long measurement

---
 test-bench-vionic.org | 76 +++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 74 insertions(+), 2 deletions(-)

diff --git a/test-bench-vionic.org b/test-bench-vionic.org
index 6332d4c..eb5769a 100644
--- a/test-bench-vionic.org
+++ b/test-bench-vionic.org
@@ -105,7 +105,7 @@ The characteristics as advertise in the manual as well as our specifications are
 | <l>               |     <c>      |       <c>       |
 | *Characteristics* |   *Manual*   | *Specification* |
 |-------------------+--------------+-----------------|
-| Time Delay        |              |    < 0.5 ms     |
+| Time Delay        |   < 10 ns    |    < 0.5 ms     |
 | Bandwidth         |  > 500 kHz   |     > 5 kHz     |
 | Noise             | < 1.6 nm rms |   < 50 nm rms   |
 | Linearity         | < +/- 15 nm  |                 |
@@ -193,11 +193,84 @@ addpath('./mat/');
 #+end_src
 
 ** TODO Thermal drifts
+
 - [ ] picture of the setup
 - [ ] long thermal drifts
+- [ ] Identification of the drifts (exponential fit)
 - [ ] once stabilize, look at the noise
 - [ ] compute low frequency ASD (may still be thermal drifts of the mechanics and not noise)
 
+#+begin_src matlab
+enc_l = load('mat/noise_meas_40h_200Hz_1.mat', 't', 'x');
+enc_l.x = enc_l.x - mean(enc_l.x(enc_l.t < 1)); % Start at zero displacement
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(enc_l.t/3600, 1e9*enc_l.x, '-');
+hold off;
+xlabel('Time [h]');
+ylabel('Displacement [nm]');
+#+end_src
+
+Exponential fit
+#+begin_src matlab
+f = @(b,x) b(1)*(1 - exp(-x/b(2)));
+
+y_cur = enc_l.x;
+t_cur = end_l.t;
+
+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
+
+Comparison of the data and exponential fit
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(enc_l.t/60/60, 1e9*enc_l.x);
+plot(enc_l.t/60/60, 1e9*f(B, enc_l.t));
+hold off;
+xlim([0, 17.5])
+xlabel('Time [h]'); ylabel('Displacement [nm]');
+#+end_src
+
+Let's get only the data once it is stabilized
+#+begin_src matlab
+x_stab = enc_l.x(enc_l.t > 20*3600);
+x_stab = x_stab - mean(x_stab);
+t_stab = enc_l.t(enc_l.t > 20*3600);
+x_stab = x_stab - x_stab(1);
+#+end_src
+
+#+begin_src matlab :exports none
+% Compute sampling Frequency
+Ts = (enc{1}.t(end) - enc{1}.t(1))/(length(enc{1}.t)-1);
+Fs = 1/Ts;
+
+% Hannning Windows
+win = hanning(ceil(60/Ts));
+
+[pxx, f] = pwelch(x_stab, win, [], [], Fs);
+#+end_src
+
+#+begin_src matlab :exports none
+figure;
+hold on;
+plot(f, sqrt(pxx))
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD [$m/\sqrt{Hz}$]');
+% xlim([10, Fs/2]);
+% ylim([1e-11, 1e-10]);
+#+end_src
+
 ** Time Domain signals
 First we load the data.
 #+begin_src matlab :exports none
@@ -302,7 +375,6 @@ exportFig('figs/vionic_noise_asd.pdf', 'width', 'wide', 'height', 'normal');
 [[file:figs/vionic_noise_asd.png]]
 
 ** Noise Model
-
 Let's create a transfer function that approximate the measured noise of the encoder.
 #+begin_src matlab
 Gn_e = 1.8e-11/(1 + s/2/pi/1e4);