From 3a8762ae9a3f0769a3224062f410e94bc6f9f38d Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Wed, 9 Sep 2020 21:13:39 +0200 Subject: [PATCH] Update orgmode and html page --- index.html | 920 ++++++++++++++++++++++++++++++++++++---- index.org | 1087 ++++++++++++++++++++++++++++++++++++++++++++---- test-bench.png | Bin 17851 -> 0 bytes test-bench.svg | 312 -------------- 4 files changed, 1853 insertions(+), 466 deletions(-) delete mode 100644 test-bench.png delete mode 100644 test-bench.svg diff --git a/index.html b/index.html index c9414b2..78b9a3e 100644 --- a/index.html +++ b/index.html @@ -3,7 +3,7 @@ "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"> - + Sensor Fusion - Test Bench @@ -15,6 +15,14 @@ + +
@@ -23,27 +31,81 @@

Table of Contents

-
-

1 Experimental Setup

+
+

1 Experimental Setup

- +

+The goal of this experimental setup is to experimentally merge inertial sensors. +

+

+To merge the sensors, optimal and robust complementary filters are designed. +

+ +

+The inertial sensors used are shown in Table +

+ +
+@@ -52,8 +114,8 @@ - - + + @@ -64,79 +126,268 @@ - +
Table 1: Inertial Sensors used
  TypeModel
GeophoneMark Product L4C - VerticalMark Product L-22 - Vertical
-
-
-
-

2 Huddle Test

-
-
-
-

2.1 Load Data

-
-
-
load('./mat/huddle_test.mat', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 't');
-dt = t(2) - t(1);
-
-
-
-
-
-

2.2 Data

-
-
-
acc_1 = acc_1 - mean(acc_1);
-acc_2 = acc_2 - mean(acc_2);
-geo_1 = geo_1 - mean(geo_1);
-geo_2 = geo_2 - mean(geo_2);
-
-
-
-
+ + + +++ ++ + + + + + + + + + + + + + + + + + + + + + + +
Table 2: Accelerometer (393B05) Specifications
SpecificationValue
Sensitivity1.02 [V/(m/s2)]
Resonant Frequency> 2.5 [kHz]
Resolution (1 to 10kHz)0.00004 [m/s2 rms]
+ + + + +++ ++ + + + + + + + + + + + + + + + + + +
Table 3: Geophone (L22) Specifications
SpecificationValue
SensitivityTo be measured [V/(m/s)]
Resonant Frequency2 [Hz]
-
-

2.3 Scale Data

-

-From raw data to estimated velocity. -This takes into account the sensibility of the sensor and possible integration to go from acceleration to velocity. +The ADC used are the IO131 Speedgoat module (link) with a 16bit resolution over +/- 10V.

-
-
G0 = 1.02; % [V/(m/s2)]
+

+The geophone signals are amplified using a DLPVA-100-B-D voltage amplified from Femto (link).
+The force sensor signal is amplified using a Low Noise Voltage Preamplifier from Ametek (link).
+

 
-G_acc = tf(G0);
-
+

+Geophone electronics: +

+
    +
  • gain: 10 (20dB)
    • +
  • low pass filter: 1.5Hz
    • +
  • hifh pass filter: 100kHz (2nd order)
    • +
+ +

+Force Sensor electronics: +

+
    +
  • gain: 10 (20dB)
    • +
  • low pass filter: 1st order at 3Hz
    • +
  • high pass filter: 1st order at 30kHz
    • +
+
-
-
T = 276;
-xi = 0.5;
-w = 2*pi;
-
-G_geo = -T*s^2/(s^2 + 2*xi*w*s + w^2);
-
+
+

2 Huddle Test

+
+

+The goal here is to measure the noise of the inertial sensors. +Is also permits to measure the motion level when nothing is actuated. +

+
+

2.1 Load Data

+
-
acc_1 = lsim(inv(G_acc), acc_1, t);
-acc_2 = lsim(inv(G_acc), acc_2, t);
-geo_1 = lsim(inv(G_geo), geo_1, t);
-geo_2 = lsim(inv(G_geo), geo_2, t);
+
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
-
-

2.4 Compare Time Domain Signals

+
+

2.2 Detrend Data

+
+
+
ht.d = detrend(ht.d, 0); % [m]
+ht.acc_1 = detrend(ht.acc_1, 0); % [V]
+ht.acc_2 = detrend(ht.acc_2, 0); % [V]
+ht.geo_1 = detrend(ht.geo_1, 0); % [V]
+ht.geo_2 = detrend(ht.geo_2, 0); % [V]
+ht.f_meas = detrend(ht.f_meas, 0); % [V]
+
+
+
+
+ +
+

2.3 Compute PSD

+
+

+We first define the parameters for the frequency domain analysis. +

+
+
Ts = t(2) - t(1); % [s]
+Fs = 1/Ts; % [Hz]
+
+win = hanning(ceil(1*Fs));
+
+
+ +

+Then we compute the Power Spectral Density using pwelch function. +

+
+
[p_d,     f] = pwelch(ht.d, win, [], [], 1/Ts);
+[p_acc1,  ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
+[p_acc2,  ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
+[p_geo1,  ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
+[p_geo2,  ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
+[p_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts);
+
+
+
+
+ +
+

2.4 Sensor Noise in Volts

+
[coh_acc, ~] = mscohere(ht.acc_1, ht.acc_2, win, [], [], Fs);
+[coh_geo, ~] = mscohere(ht.geo_1, ht.geo_2, win, [], [], Fs);
+
+
+ +
+
pN_acc = p_acc1.*(1 - coh_acc);
+pN_geo = p_geo1.*(1 - coh_geo);
+
+
+ +

+PSD of the ADC quantization noise. +

+
+
Sq = (20/2^16)^2/(12*Fs);
+
+
+ +
+
figure;
+hold on;
+plot(f, sqrt(pN_acc), '-', 'DisplayName', 'Accelerometers');
+plot(f, sqrt(pN_geo), '-', 'DisplayName', 'Geophones');
+plot(f, ones(size(f))*sqrt(Sq), '-', 'DisplayName', 'ADC');
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the Measurement Noise $[V/\sqrt{Hz}]$');
+xlim([1, 5000]);
+legend('location', 'northeast');
+
+
+
+
+
+ +
+

3 After identification

+
+
+
+

3.1 Scale Data

+
+

+Let’s use a model of the accelerometer and geophone to compute the motion from the measured voltage. +

+
+
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
+
+
+ +
+
G_geo = 120*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [V/(m/s)]
+
+
+ +
+
figure;
+hold on;
+set(gca, 'ColorOrderIndex', 1);
+plot(f, sqrt(p_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
+     'DisplayName', 'Accelerometer');
+set(gca, 'ColorOrderIndex', 1);
+plot(f, sqrt(p_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
+     'HandleVisibility', 'off');
+set(gca, 'ColorOrderIndex', 2);
+plot(f, sqrt(p_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
+     'DisplayName', 'Geophone');
+set(gca, 'ColorOrderIndex', 2);
+plot(f, sqrt(p_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
+     'HandleVisibility', 'off');
+set(gca, 'ColorOrderIndex', 3);
+plot(f, sqrt(p_d), 'DisplayName', 'Interferometer');
+hold off;
+set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
+ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
+title('Huddle Test')
+legend();
+
+
+
+
+ +
+

3.2 sdlkfj

+
+
+
acc_1 = lsim(inv(G_acc), ht.acc_1, ht.t);
+acc_2 = lsim(inv(G_acc), ht.acc_2, ht.t);
+
+geo_1 = lsim(inv(G_geo), ht.geo_1, ht.t);
+geo_2 = lsim(inv(G_geo), ht.geo_2, ht.t);
+
+
+
+
+ +
+

3.3 Compare Time Domain Signals

+
+
figure;
 hold on;
 plot(t, acc_1);
@@ -149,9 +400,9 @@ hold off;
-
-

2.5 Compute PSD

-
+
+

3.4 Compute PSD

+

We first define the parameters for the frequency domain analysis.

@@ -201,9 +452,9 @@ xlim([1, 5000]);
-
-

2.6 Dynamical Uncertainty

-
+
+

3.5 Dynamical Uncertainty

+
[T_acc, ~] = tfestimate(acc_1, acc_2, win, [], [], Fs);
 [T_geo, ~] = tfestimate(geo_1, geo_2, win, [], [], Fs);
@@ -212,9 +463,37 @@ xlim([1, 5000]);
-
-

2.7 Sensor Noise

-
+
+

3.6 ADC Noise

+
+

+Let’s note: +

+
    +
  • \(\Delta V\) the ADC range in [V]
    • +
  • \(n\) the number of bits
    • +
  • \(q = \frac{\Delta V}{2^n}\) the quantization in [V]
    • +
  • \(f_N\) the sampling frequency in [Hz]
    • +
+ +

+The Power Spectral Density of the quantization noise is then: +

+\begin{equation} + S_Q = \frac{q^2}{12 f_N} \quad [V^2/Hz] +\end{equation} + +
+
Fs = 1/dt;
+Sq = (20/2^16)^2/(12*Fs);
+
+
+
+
+ +
+

3.7 Sensor Noise

+
[coh_acc, ~] = mscohere(acc_1, acc_2, win, [], [], Fs);
 [coh_geo, ~] = mscohere(geo_1, geo_2, win, [], [], Fs);
@@ -242,10 +521,499 @@ legend('location', 'northeast');
+ +
+

4 Sensor Dynamics

+
+

+Thanks to the interferometer, it is possible to compute the transfer function from the mass displacement to the voltage generated by the inertial sensors. +This permits to estimate the sensor dynamics and to calibrate the sensors. +

+
+ +
+

4.1 Load Data

+
+
+
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+id_ol = load('./mat/identification_noise_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+
+
+
+
+ +
+

4.2 Time Domain Signals

+
+

+Excitation signal: noise. +

+ + +
+

first_exc_signal_time.png +

+

Figure 1: Excitation signal used for the first identification

+
+
+
+ +
+

4.3 Identification of the IFF Plant

+
+
+
+

4.3.1 Experimental Data

+
+
+
Ts = ht.t(2) - ht.t(1);
+win = hann(ceil(10/Ts));
+
+
+ +
+
[tf_fmeas_est, f] = tfestimate(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts); % [V/m]
+[co_fmeas_est, ~] = mscohere(  id_ol.u, id_ol.f_meas, win, [], [], 1/Ts);
+
+
+ + +
+

iff_identification_coh.png +

+

Figure 2: Coherence for the identification of the IFF plant

+
+ + +
+

iff_identification_bode_plot.png +

+

Figure 3: Bode plot of the identified IFF plant

+
+
+
+ +
+

4.3.2 Model of the IFF Plant

+
+
+
wz = 2*pi*102;
+xi_z = 0.01;
+wp = 2*pi*239.4;
+xi_p = 0.015;
+
+Giff = 2.2*(s^2 + 2*xi_z*s*wz + wz^2)/(s^2 + 2*xi_p*s*wp + wp^2) * ... % Dynamics
+       10*(s/3/pi/(1 + s/3/pi)) * ... % Low pass filter and gain of the voltage amplifier
+       exp(-Ts*s); % Time delay induced by ADC/DAC
+
+
+ + +
+

iff_plant_model.png +

+

Figure 4: IFF Plant + Model

+
+
+
+ +
+

4.3.3 Root Locus and optimal Controller

+
+ +
+

iff_root_locus.png +

+

Figure 5: Root Locus for the IFF control

+
+ +

+The controller that yield maximum damping is: +

+
+
Kiff_opt = 102/(s + 2*pi*2);
+
+
+
+
+
+ +
+

4.4 Identification of Sensor Dynamics with IFF activated

+
+
+
+

4.4.1 Signals

+
+

+A new identification is performed with the resonance damped. +It helps to not induce too much motion at the resonance and damage the actuator. +

+
+
id_cl = load('./mat/identification_noise_iff_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+
+
+
+
+ +
+

4.4.2 Verification of the achievable damping

+
+
+
[tf_G_ol_est, f] = tfestimate(id_ol.u, id_ol.d, win, [], [], 1/Ts);
+[co_G_ol_est, ~] = mscohere(  id_ol.u, id_ol.d, win, [], [], 1/Ts);
+[tf_G_cl_est, ~] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts);
+[co_G_cl_est, ~] = mscohere(  id_cl.u, id_cl.d, win, [], [], 1/Ts);
+
+
+ + +
+

Gd_identification_iff_coherence.png +

+

Figure 6: Coherence for the transfer function from F to d, with and without IFF

+
+ +

+Don’t really understand the low frequency behavior. +

+ +
+

Gd_identification_iff_bode_plot.png +

+

Figure 7: Coherence for the transfer function from F to d, with and without IFF

+
+
+
+
+ +
+

4.5 Generate the excitation signal

+
+
+
+

4.5.1 Requirements

+
+

+The requirements on the excitation signal is: +

+
    +
  • General much larger motion that the measured motion during the huddle test
    • +
  • Don’t damage the actuator
    • +
+ +

+To determine the perfect voltage signal to be generated, we need two things: +

+
    +
  • the transfer function from voltage to mass displacement
    • +
  • the PSD of the measured motion by the inertial sensors
    • +
  • not saturate the sensor signals
    • +
  • provide enough signal/noise ratio (good coherence) in the frequency band of interest (~0.5Hz to 3kHz)
    • +
+
+
+ +
+

4.5.2 Transfer function from excitation signal to displacement

+
+

+Let’s first estimate the transfer function from the excitation signal in [V] to the generated displacement in [m] as measured by the inteferometer. +

+ +
+
id_cl = load('./mat/identification_noise_iff_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+
+
+ +
+
win = hann(ceil(10/Ts));
+
+[tf_G_cl_est, f] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts);
+[co_G_cl_est, ~] = mscohere(  id_cl.u, id_cl.d, win, [], [], 1/Ts);
+
+
+ +

+Approximate transfer function from voltage output to generated displacement when IFF is used, in [m/V]. +

+
+
G_d_est = -5e-6*(2*pi*230)^2/(s^2 + 2*0.3*2*pi*240*s + (2*pi*240)^2);
+
+
+ + +
+

Gd_plant_estimation.png +

+

Figure 8: Estimation of the transfer function from the excitation signal to the generated displacement

+
+
+
+ +
+

4.5.3 Motion measured during Huddle test

+
+

+We now compute the PSD of the measured motion by the inertial sensors during the huddle test. +

+
+
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+
+
+ +
+
[p_d, f] = pwelch(ht.d, win, [], [], 1/Ts);
+[p_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
+[p_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
+[p_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
+[p_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
+
+
+ +

+Using an estimated model of the sensor dynamics from the documentation of the sensors, we can compute the ASD of the motion in \(m/\sqrt{Hz}\) measured by the sensors. +

+
+
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
+G_geo = -120*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [V/(m/s)]
+
+
+ + +
+

huddle_test_psd_motion.png +

+

Figure 9: ASD of the motion measured by the sensors

+
+ +

+From the ASD of the motion measured by the sensors, we can create an excitation signal that will generate much motion motion that the motion under no excitation. +

+ +

+We create G_exc that corresponds to the wanted generated motion. +

+
+
G_exc = 0.2e-6/(1 + s/2/pi/2)/(1 + s/2/pi/50);
+
+
+ +

+And we create a time domain signal y_d that have the spectral density described by G_exc. +

+
+
Fs = 1/Ts;
+t = 0:Ts:180; % Time Vector [s]
+u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one
+
+y_d = lsim(G_exc, u, t);
+
+
+ +
+
[pxx, ~] = pwelch(y_d, win, 0, [], Fs);
+
+
+ + +
+

comp_huddle_test_excited_motion_psd.png +

+

Figure 10: Comparison of the ASD of the motion during Huddle and the wanted generated motion

+
+ + +

+We can now generate the voltage signal that will generate the wanted motion. +

+
+
y_v = lsim(G_exc * ... % from unit PSD to shaped PSD
+           (1 + s/2/pi/50) * ... % Inverse of pre-filter included in the Simulink file
+           1/G_d_est * ... % Wanted displacement => required voltage
+           1/(1 + s/2/pi/5e3), ... %  Add some high frequency filtering
+           u, t);
+
+
+ + +
+

optimal_exc_signal_time.png +

+

Figure 11: Generated excitation signal

+
+
+
+
+ +
+

4.6 Identification of the Inertial Sensors Dynamics

+
+
+
+

4.6.1 Load Data

+
+
+
id = load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
+
+
+ +
+
ht.d = detrend(ht.d, 0);
+ht.acc_1 = detrend(ht.acc_1, 0);
+ht.acc_2 = detrend(ht.acc_2, 0);
+ht.geo_1 = detrend(ht.geo_1, 0);
+ht.geo_2 = detrend(ht.geo_2, 0);
+ht.f_meas = detrend(ht.f_meas, 0);
+
+
+ +
+
id.d = detrend(id.d, 0);
+id.acc_1 = detrend(id.acc_1, 0);
+id.acc_2 = detrend(id.acc_2, 0);
+id.geo_1 = detrend(id.geo_1, 0);
+id.geo_2 = detrend(id.geo_2, 0);
+id.f_meas = detrend(id.f_meas, 0);
+
+
+
+
+ +
+

4.6.2 Compare PSD during Huddle and and during identification

+
+
+
Ts = ht.t(2) - ht.t(1);
+win = hann(ceil(10/Ts));
+
+
+ +
+
[p_id_d, f] = pwelch(id.d, win, [], [], 1/Ts);
+[p_id_acc1, ~] = pwelch(id.acc_1, win, [], [], 1/Ts);
+[p_id_acc2, ~] = pwelch(id.acc_2, win, [], [], 1/Ts);
+[p_id_geo1, ~] = pwelch(id.geo_1, win, [], [], 1/Ts);
+[p_id_geo2, ~] = pwelch(id.geo_2, win, [], [], 1/Ts);
+[p_id_fmeas, ~] = pwelch(id.f_meas, win, [], [], 1/Ts);
+
+
+ +
+
[p_ht_d, ~] = pwelch(ht.d, win, [], [], 1/Ts);
+[p_ht_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
+[p_ht_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
+[p_ht_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
+[p_ht_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
+[p_ht_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts);
+
+
+ + +
+

comp_psd_huddle_test_identification.png +

+

Figure 12: Comparison of the PSD of the measured motion during the Huddle test and during the identification

+
+
+
+ +
+

4.6.3 Compute transfer functions

+
+
+
[tf_acc1_est, f] = tfestimate(id.d, id.acc_1, win, [], [], 1/Ts);
+[co_acc1_est, ~] = mscohere(  id.d, id.acc_1, win, [], [], 1/Ts);
+[tf_acc2_est, ~] = tfestimate(id.d, id.acc_2, win, [], [], 1/Ts);
+[co_acc2_est, ~] = mscohere(  id.d, id.acc_2, win, [], [], 1/Ts);
+
+[tf_geo1_est, ~] = tfestimate(id.d, id.geo_1, win, [], [], 1/Ts);
+[co_geo1_est, ~] = mscohere(  id.d, id.geo_1, win, [], [], 1/Ts);
+[tf_geo2_est, ~] = tfestimate(id.d, id.geo_2, win, [], [], 1/Ts);
+[co_geo2_est, ~] = mscohere(  id.d, id.geo_2, win, [], [], 1/Ts);
+
+
+ + +
+

id_sensor_dynamics_coherence.png +

+

Figure 13: Coherence for the estimation of the sensor dynamics

+
+ +

+Model of the inertial sensors: +

+
+
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
+G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)]
+
+
+ + +
+

id_sensor_dynamics_accelerometers.png +

+

Figure 14: Identified dynamics of the accelerometers

+
+ + + +
+

id_sensor_dynamics_geophones.png +

+

Figure 15: Identified dynamics of the geophones

+
+
+
+
+ +
+

4.7 Compare Time domain Estimation of the displacement

+
+

+Let’s compare the measured accelerations instead of displacement (no integration). +

+ +
+
G_lpf = 1/(1 + s/2/pi/5e3);
+
+acc1_a = lsim(1/G_acc*G_lpf, id.acc_1, id.t);
+acc2_a = lsim(1/G_acc*G_lpf, id.acc_2, id.t);
+
+
+ +
+
geo1_a = lsim(1/G_geo*s*G_lpf, id.geo_1, id.t);
+geo2_a = lsim(1/G_geo*s*G_lpf, id.geo_2, id.t);
+
+
+ +
+
int_a = lsim(s^2*G_lpf*G_lpf, id.d, id.t);
+
+
+ +
+
figure;
+hold on;
+plot(id.t, int_a);
+plot(id.t, acc1_a);
+plot(id.t, acc2_a);
+plot(id.t, geo1_a);
+plot(id.t, geo2_a);
+hold off;
+xlabel('Time [s]'); ylabel('Acceleration [m]');
+
+
+
+
+

Author: Dehaeze Thomas

-

Created: 2020-08-31 lun. 16:09

+

Created: 2020-09-09 mer. 20:38

diff --git a/index.org b/index.org index 7b86462..d7dbd5d 100644 --- a/index.org +++ b/index.org @@ -36,13 +36,56 @@ * Experimental Setup -| | | +The goal of this experimental setup is to experimentally merge inertial sensors. + +To merge the sensors, optimal and robust complementary filters are designed. + +The inertial sensors used are shown in Table + +#+name: tab:inertial_sensors +#+caption: Inertial Sensors used +| *Type* | *Model* | |---------------+------------------------------| | Accelerometer | PCB 393B05 - Vertical ([[https://www.pcb.com/products?m=393B05][link]]) | -| Geophone | Mark Product L4C - Vertical | +| Geophone | Mark Product L-22 - Vertical | + + +#+name: tab:393B05_spec +#+caption: Accelerometer (393B05) Specifications +| *Specification* | *Value* | +|-------------------------+--------------------| +| Sensitivity | 1.02 [V/(m/s2)] | +| Resonant Frequency | > 2.5 [kHz] | +| Resolution (1 to 10kHz) | 0.00004 [m/s2 rms] | + +#+name: tab:L22_spec +#+caption: Geophone (L22) Specifications +| *Specification* | *Value* | +|-------------------------+--------------------------| +| Sensitivity | To be measured [V/(m/s)] | +| Resonant Frequency | 2 [Hz] | + +The ADC used are the IO131 Speedgoat module ([[https://www.speedgoat.com/products/io-connectivity-analog-io131][link]]) with a 16bit resolution over +/- 10V. + +The geophone signals are amplified using a DLPVA-100-B-D voltage amplified from Femto ([[https://www.femto.de/en/products/voltage-amplifiers/variable-gain-100-khz-dlpva.html][link]]). +The force sensor signal is amplified using a Low Noise Voltage Preamplifier from Ametek ([[https://www.ameteksi.com/support-center/legacy-products/signal-recovery-legacy/5113-low-noise-preamplifier][link]]). + +Geophone electronics: +- gain: 10 (20dB) +- low pass filter: 1.5Hz +- hifh pass filter: 100kHz (2nd order) + +Force Sensor electronics: +- gain: 10 (20dB) +- low pass filter: 1st order at 3Hz +- high pass filter: 1st order at 30kHz * Huddle Test -** Matlab Init :noexport:ignore: +** Introduction :ignore: +The goal here is to measure the noise of the inertial sensors. +Is also permits to measure the motion level when nothing is actuated. + +** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) <> #+end_src @@ -53,99 +96,979 @@ ** Load Data #+begin_src matlab - load('./mat/huddle_test.mat', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 't'); - dt = t(2) - t(1); + ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); #+end_src -** Data +** Detrend Data #+begin_src matlab - acc_1 = acc_1 - mean(acc_1); - acc_2 = acc_2 - mean(acc_2); - geo_1 = geo_1 - mean(geo_1); - geo_2 = geo_2 - mean(geo_2); -#+end_src - -** Scale Data -From raw data to estimated velocity. -This takes into account the sensibility of the sensor and possible integration to go from acceleration to velocity. - -#+begin_src matlab - G0 = 1.02; % [V/(m/s2)] - - G_acc = tf(G0); -#+end_src - -#+begin_src matlab - T = 276; - xi = 0.5; - w = 2*pi; - - G_geo = -T*s^2/(s^2 + 2*xi*w*s + w^2); -#+end_src - -#+begin_src matlab - acc_1 = lsim(inv(G_acc), acc_1, t); - acc_2 = lsim(inv(G_acc), acc_2, t); - geo_1 = lsim(inv(G_geo), geo_1, t); - geo_2 = lsim(inv(G_geo), geo_2, t); -#+end_src - -** Compare Time Domain Signals -#+begin_src matlab - figure; - hold on; - plot(t, acc_1); - plot(t, acc_2); - plot(t, geo_1); - plot(t, geo_2); - hold off; + ht.d = detrend(ht.d, 0); % [m] + ht.acc_1 = detrend(ht.acc_1, 0); % [V] + ht.acc_2 = detrend(ht.acc_2, 0); % [V] + ht.geo_1 = detrend(ht.geo_1, 0); % [V] + ht.geo_2 = detrend(ht.geo_2, 0); % [V] + ht.f_meas = detrend(ht.f_meas, 0); % [V] #+end_src ** Compute PSD We first define the parameters for the frequency domain analysis. -#+begin_src matlab :results none - Fs = 1/dt; % [Hz] +#+begin_src matlab + Ts = ht.t(2) - ht.t(1); % [s] + Fs = 1/Ts; % [Hz] win = hanning(ceil(1*Fs)); #+end_src Then we compute the Power Spectral Density using =pwelch= function. -#+begin_src matlab :results none - [p_acc_1, f] = pwelch(acc_1, win, [], [], Fs); - [p_acc_2, ~] = pwelch(acc_2, win, [], [], Fs); - [p_geo_1, ~] = pwelch(geo_1, win, [], [], Fs); - [p_geo_2, ~] = pwelch(geo_2, win, [], [], Fs); +#+begin_src matlab + [p_d, f] = pwelch(ht.d, win, [], [], 1/Ts); + [p_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts); + [p_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts); + [p_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts); + [p_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts); + [p_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(1, 2, 1); + hold on; + plot(f, sqrt(p_acc1)); + plot(f, sqrt(p_acc2)); + hold off; + set(gca, 'xscale', 'log'); + set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('ASD Accelerometers $[V/\sqrt{Hz}]$') + xlim([1, 5000]); + + ax2 = subplot(1, 2, 2); + hold on; + plot(f, sqrt(p_geo1)); + plot(f, sqrt(p_geo2)); + hold off; + set(gca, 'xscale', 'log'); + set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('ASD Geophones $[V/\sqrt{Hz}]$') + xlim([1, 5000]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace +exportFig('figs/huddle_test_measured_voltage_inertial_sensors.pdf', 'width', 'full', 'height', 'tall'); +#+end_src + +#+name: fig:huddle_test_measured_voltage_inertial_sensors +#+caption: ASD of the measured voltage from the inertial sensors during the Huddle test +#+RESULTS: +[[file:figs/huddle_test_measured_voltage_inertial_sensors.png]] + +** Sensor Noise in Volts +#+begin_src matlab + [coh_acc, ~] = mscohere(ht.acc_1, ht.acc_2, win, [], [], Fs); + [coh_geo, ~] = mscohere(ht.geo_1, ht.geo_2, win, [], [], Fs); +#+end_src + +#+begin_src matlab + pN_acc = p_acc1.*(1 - coh_acc); + pN_geo = p_geo1.*(1 - coh_geo); +#+end_src + +PSD of the ADC quantization noise. +#+begin_src matlab + Sq = (20/2^16)^2/(12*Fs); #+end_src #+begin_src matlab :results none figure; hold on; - plot(f, sqrt(p_acc_1)); - plot(f, sqrt(p_acc_2)); + plot(f, sqrt(pN_acc), '-', 'DisplayName', 'Accelerometers'); + plot(f, sqrt(pN_geo), '-', 'DisplayName', 'Geophones'); + plot(f, ones(size(f))*sqrt(Sq), '-', 'DisplayName', 'ADC'); hold off; - set(gca, 'xscale', 'log'); - set(gca, 'yscale', 'log'); - xlabel('Frequency [Hz]'); ylabel('ASD Accelerometers $\left[\frac{m/s}{\sqrt{Hz}}\right]$') + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('ASD of the Measurement Noise $[V/\sqrt{Hz}]$'); xlim([1, 5000]); + legend('location', 'northeast'); +#+end_src + +* Sensor Dynamics +** Introduction :ignore: +Thanks to the interferometer, it is possible to compute the transfer function from the mass displacement to the voltage generated by the inertial sensors. +This permits to estimate the sensor dynamics and to calibrate the sensors. + +** Matlab Init :noexport:ignore: +#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) + <> +#+end_src + +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +** Load Data +#+begin_src matlab + ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); + id_ol = load('./mat/identification_noise_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab :exports none + ht.d = detrend(ht.d, 0); % [m] + ht.acc_1 = detrend(ht.acc_1, 0); % [V] + ht.acc_2 = detrend(ht.acc_2, 0); % [V] + ht.geo_1 = detrend(ht.geo_1, 0); % [V] + ht.geo_2 = detrend(ht.geo_2, 0); % [V] + ht.f_meas = detrend(ht.f_meas, 0); % [V] +#+end_src + +#+begin_src matlab :exports none + id_ol.d = detrend(id_ol.d, 0); + id_ol.acc_1 = detrend(id_ol.acc_1, 0); + id_ol.acc_2 = detrend(id_ol.acc_2, 0); + id_ol.geo_1 = detrend(id_ol.geo_1, 0); + id_ol.geo_2 = detrend(id_ol.geo_2, 0); + id_ol.f_meas = detrend(id_ol.f_meas, 0); + id_ol.u = detrend(id_ol.u, 0); +#+end_src + +** Time Domain Signals +Excitation signal: noise. + +#+begin_src matlab :exports none + figure; + plot(id_ol.t, id_ol.u); + xlim([0 0.1]); + xlabel('Time [s]'); ylabel('Excitation Signal [V]'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/first_exc_signal_time.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:first_exc_signal_time +#+caption: Excitation signal used for the first identification +#+RESULTS: +[[file:figs/first_exc_signal_time.png]] + +** Identification of the IFF Plant +*** Experimental Data +#+begin_src matlab + Ts = ht.t(2) - ht.t(1); + win = hann(ceil(10/Ts)); +#+end_src + +#+begin_src matlab + [tf_fmeas_est, f] = tfestimate(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts); % [V/m] + [co_fmeas_est, ~] = mscohere( id_ol.u, id_ol.f_meas, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + plot(f, co_fmeas_est, '-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Coherence'); xlabel('Frequency [Hz]'); + hold off; + xlim([1, 1e3]); ylim([0, 1]) +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/iff_identification_coh.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:iff_identification_coh +#+caption: Coherence for the identification of the IFF plant +#+RESULTS: +[[file:figs/iff_identification_coh.png]] + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_fmeas_est), '-') + 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_fmeas_est), '-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2], 'x'); + xlim([1, 1e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/iff_identification_bode_plot.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:iff_identification_bode_plot +#+caption: Bode plot of the identified IFF plant +#+RESULTS: +[[file:figs/iff_identification_bode_plot.png]] + +*** Model of the IFF Plant +#+begin_src matlab + wz = 2*pi*102; + xi_z = 0.01; + wp = 2*pi*239.4; + xi_p = 0.015; + + Giff = 2.2*(s^2 + 2*xi_z*s*wz + wz^2)/(s^2 + 2*xi_p*s*wp + wp^2) * ... % Dynamics + 10*(s/3/pi/(1 + s/3/pi)) * ... % Low pass filter and gain of the voltage amplifier + exp(-Ts*s); % Time delay induced by ADC/DAC +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_fmeas_est), '.') + plot(f, abs(squeeze(freqresp(Giff, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude [V/V]'); set(gca, 'XTickLabel',[]); + hold off; + ylim([1e-2, 1e3]) + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_fmeas_est), '.') + plot(f, 180/pi*angle(squeeze(freqresp(Giff, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + + linkaxes([ax1,ax2], 'x'); + xlim([0.5, 5e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/iff_plant_model.pdf', 'width', 'large', 'height', 'full'); +#+end_src + +#+name: fig:iff_plant_model +#+caption: IFF Plant + Model +#+RESULTS: +[[file:figs/iff_plant_model.png]] + +*** Root Locus and optimal Controller + +#+begin_src matlab :exports none + gains = logspace(0, 5, 1000); + + figure; + hold on; + plot(real(pole(Giff)), imag(pole(Giff)), 'kx'); + plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko'); + for i = 1:length(gains) + cl_poles = pole(feedback(Giff, gains(i)/(s + 2*pi*2))); + plot(real(cl_poles), imag(cl_poles), 'k.'); + end + cl_poles = pole(feedback(Giff, 102/(s + 2*pi*2))); + plot(real(cl_poles), imag(cl_poles), 'rx'); + ylim([0, 1800]); + xlim([-1600,200]); + xlabel('Real Part') + ylabel('Imaginary Part') + axis square +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/iff_root_locus.pdf', 'width', 'wide', 'height', 'tall'); +#+end_src + +#+name: fig:iff_root_locus +#+caption: Root Locus for the IFF control +#+RESULTS: +[[file:figs/iff_root_locus.png]] + +The controller that yield maximum damping is: +#+begin_src matlab + Kiff_opt = 102/(s + 2*pi*2); +#+end_src + +** Identification of Sensor Dynamics with IFF activated +*** Signals +A new identification is performed with the resonance damped. +It helps to not induce too much motion at the resonance and damage the actuator. +#+begin_src matlab + id_cl = load('./mat/identification_noise_iff_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab :exports none + id_cl.d = detrend(id_cl.d, 0); + id_cl.acc_1 = detrend(id_cl.acc_1, 0); + id_cl.acc_2 = detrend(id_cl.acc_2, 0); + id_cl.geo_1 = detrend(id_cl.geo_1, 0); + id_cl.geo_2 = detrend(id_cl.geo_2, 0); + id_cl.f_meas = detrend(id_cl.f_meas, 0); + id_cl.u = detrend(id_cl.u, 0); +#+end_src + +*** Verification of the achievable damping +#+begin_src matlab + [tf_G_ol_est, f] = tfestimate(id_ol.u, id_ol.d, win, [], [], 1/Ts); + [co_G_ol_est, ~] = mscohere( id_ol.u, id_ol.d, win, [], [], 1/Ts); + [tf_G_cl_est, ~] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts); + [co_G_cl_est, ~] = mscohere( id_cl.u, id_cl.d, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + plot(f, co_G_ol_est, '-', 'DisplayName', 'OL') + plot(f, co_G_cl_est, '-', 'DisplayName', 'IFF') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Coherence'); xlabel('Frequency [Hz]'); + hold off; + xlim([1, 1e3]); ylim([0, 1]) + legend('location', 'southwest'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/Gd_identification_iff_coherence.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:Gd_identification_iff_coherence +#+caption: Coherence for the transfer function from F to d, with and without IFF +#+RESULTS: +[[file:figs/Gd_identification_iff_coherence.png]] + +Don't really understand the low frequency behavior. +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_G_ol_est), '-', 'DisplayName', 'OL') + plot(f, abs(tf_G_cl_est), '-', 'DisplayName', 'IFF') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude [m/V]'); set(gca, 'XTickLabel',[]); + hold off; + legend('location', 'northeast'); + ylim([2e-7, 2e-4]); + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_G_ol_est), '-') + plot(f, 180/pi*angle(tf_G_cl_est), '-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + + linkaxes([ax1,ax2], 'x'); + xlim([1, 1e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/Gd_identification_iff_bode_plot.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:Gd_identification_iff_bode_plot +#+caption: Coherence for the transfer function from F to d, with and without IFF +#+RESULTS: +[[file:figs/Gd_identification_iff_bode_plot.png]] + +** Generate the excitation signal +*** Requirements +The requirements on the excitation signal is: +- General much larger motion that the measured motion during the huddle test +- Don't damage the actuator + +To determine the perfect voltage signal to be generated, we need two things: +- the transfer function from voltage to mass displacement +- the PSD of the measured motion by the inertial sensors +- not saturate the sensor signals +- provide enough signal/noise ratio (good coherence) in the frequency band of interest (~0.5Hz to 3kHz) + +*** Transfer function from excitation signal to displacement +Let's first estimate the transfer function from the excitation signal in [V] to the generated displacement in [m] as measured by the inteferometer. + +#+begin_src matlab + id_cl = load('./mat/identification_noise_iff_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab :exports none + id_cl.d = detrend(id_cl.d, 0); + id_cl.acc_1 = detrend(id_cl.acc_1, 0); + id_cl.acc_2 = detrend(id_cl.acc_2, 0); + id_cl.geo_1 = detrend(id_cl.geo_1, 0); + id_cl.geo_2 = detrend(id_cl.geo_2, 0); + id_cl.f_meas = detrend(id_cl.f_meas, 0); + id_cl.u = detrend(id_cl.u, 0); +#+end_src + +#+begin_src matlab + win = hann(ceil(10/Ts)); + + [tf_G_cl_est, f] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts); + [co_G_cl_est, ~] = mscohere( id_cl.u, id_cl.d, win, [], [], 1/Ts); +#+end_src + +Approximate transfer function from voltage output to generated displacement when IFF is used, in [m/V]. +#+begin_src matlab + G_d_est = -5e-6*(2*pi*230)^2/(s^2 + 2*0.3*2*pi*240*s + (2*pi*240)^2); +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_G_cl_est), '-') + plot(f, abs(squeeze(freqresp(G_d_est, f, 'Hz'))), '--') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude [m/V]'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_G_cl_est), '-') + plot(f, 180/pi*angle(squeeze(freqresp(G_d_est, f, 'Hz'))), '--') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2], 'x'); + xlim([10, 1000]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/Gd_plant_estimation.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:Gd_plant_estimation +#+caption: Estimation of the transfer function from the excitation signal to the generated displacement +#+RESULTS: +[[file:figs/Gd_plant_estimation.png]] + +*** Motion measured during Huddle test +We now compute the PSD of the measured motion by the inertial sensors during the huddle test. +#+begin_src matlab + ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab :exports none + ht.d = detrend(ht.d, 0); + ht.acc_1 = detrend(ht.acc_1, 0); + ht.acc_2 = detrend(ht.acc_2, 0); + ht.geo_1 = detrend(ht.geo_1, 0); + ht.geo_2 = detrend(ht.geo_2, 0); +#+end_src + +#+begin_src matlab + [p_d, f] = pwelch(ht.d, win, [], [], 1/Ts); + [p_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts); + [p_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts); + [p_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts); + [p_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts); +#+end_src + +Using an estimated model of the sensor dynamics from the documentation of the sensors, we can compute the ASD of the motion in $m/\sqrt{Hz}$ measured by the sensors. +#+begin_src matlab + G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)] + G_geo = -120*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [V/(m/s)] +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'DisplayName', 'Accelerometer'); + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'DisplayName', 'Geophone'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 3); + plot(f, sqrt(pxx), 'k-', ... + 'DisplayName', 'Excitation'); + plot(f, sqrt(p_d), 'DisplayName', 'Interferometer'); + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]'); + title('Huddle Test') + legend(); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/huddle_test_psd_motion.pdf', 'width', 'wide', 'height', 'tall'); +#+end_src + +#+name: fig:huddle_test_psd_motion +#+caption: ASD of the motion measured by the sensors +#+RESULTS: +[[file:figs/huddle_test_psd_motion.png]] + +From the ASD of the motion measured by the sensors, we can create an excitation signal that will generate much motion motion that the motion under no excitation. + +We create =G_exc= that corresponds to the wanted generated motion. +#+begin_src matlab + G_exc = 0.2e-6/(1 + s/2/pi/2)/(1 + s/2/pi/50); +#+end_src + +And we create a time domain signal =y_d= that have the spectral density described by =G_exc=. +#+begin_src matlab + Fs = 1/Ts; + t = 0:Ts:180; % Time Vector [s] + u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one + + y_d = lsim(G_exc, u, t); +#+end_src + +#+begin_src matlab + [pxx, ~] = pwelch(y_d, win, 0, [], Fs); +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'DisplayName', 'Accelerometer'); + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'DisplayName', 'Geophone'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 3); + plot(f, sqrt(pxx), 'k-', ... + 'DisplayName', 'Excitation'); + plot(f, sqrt(p_d), 'DisplayName', 'Interferometer'); + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]'); + title('Huddle Test') + legend(); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/comp_huddle_test_excited_motion_psd.pdf', 'width', 'wide', 'height', 'tall'); +#+end_src + +#+name: fig:comp_huddle_test_excited_motion_psd +#+caption: Comparison of the ASD of the motion during Huddle and the wanted generated motion +#+RESULTS: +[[file:figs/comp_huddle_test_excited_motion_psd.png]] + + +We can now generate the voltage signal that will generate the wanted motion. +#+begin_src matlab + y_v = lsim(G_exc * ... % from unit PSD to shaped PSD + (1 + s/2/pi/50) * ... % Inverse of pre-filter included in the Simulink file + 1/G_d_est * ... % Wanted displacement => required voltage + 1/(1 + s/2/pi/5e3), ... % Add some high frequency filtering + u, t); +#+end_src + +#+begin_src matlab :exports none + figure; + plot(t, y_v) + xlabel('Time [s]'); ylabel('Voltage [V]'); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/optimal_exc_signal_time.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:optimal_exc_signal_time +#+caption: Generated excitation signal +#+RESULTS: +[[file:figs/optimal_exc_signal_time.png]] + +** Identification of the Inertial Sensors Dynamics +*** Load Data +#+begin_src matlab + id = load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); + ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab + ht.d = detrend(ht.d, 0); + ht.acc_1 = detrend(ht.acc_1, 0); + ht.acc_2 = detrend(ht.acc_2, 0); + ht.geo_1 = detrend(ht.geo_1, 0); + ht.geo_2 = detrend(ht.geo_2, 0); + ht.f_meas = detrend(ht.f_meas, 0); +#+end_src + +#+begin_src matlab + id.d = detrend(id.d, 0); + id.acc_1 = detrend(id.acc_1, 0); + id.acc_2 = detrend(id.acc_2, 0); + id.geo_1 = detrend(id.geo_1, 0); + id.geo_2 = detrend(id.geo_2, 0); + id.f_meas = detrend(id.f_meas, 0); +#+end_src + +*** Compare PSD during Huddle and and during identification +#+begin_src matlab + Ts = ht.t(2) - ht.t(1); + win = hann(ceil(10/Ts)); +#+end_src + +#+begin_src matlab + [p_id_d, f] = pwelch(id.d, win, [], [], 1/Ts); + [p_id_acc1, ~] = pwelch(id.acc_1, win, [], [], 1/Ts); + [p_id_acc2, ~] = pwelch(id.acc_2, win, [], [], 1/Ts); + [p_id_geo1, ~] = pwelch(id.geo_1, win, [], [], 1/Ts); + [p_id_geo2, ~] = pwelch(id.geo_2, win, [], [], 1/Ts); + [p_id_fmeas, ~] = pwelch(id.f_meas, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab + [p_ht_d, ~] = pwelch(ht.d, win, [], [], 1/Ts); + [p_ht_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts); + [p_ht_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts); + [p_ht_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts); + [p_ht_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts); + [p_ht_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(1, 2, 1); + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, p_ht_acc1, 'DisplayName', 'Huddle Test'); + set(gca, 'ColorOrderIndex', 1); + plot(f, p_ht_acc2, 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 2); + plot(f, p_id_acc1, 'DisplayName', 'Identification Test'); + set(gca, 'ColorOrderIndex', 2); + plot(f, p_id_acc2, 'HandleVisibility', 'off'); + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]'); + title('Huddle Test - Accelerometers') + legend(); + + ax1 = subplot(1, 2, 2); + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, p_ht_geo1, 'DisplayName', 'Huddle Test'); + set(gca, 'ColorOrderIndex', 1); + plot(f, p_ht_geo2, 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 2); + plot(f, p_id_geo1, 'DisplayName', 'Identification Test'); + set(gca, 'ColorOrderIndex', 2); + plot(f, p_id_geo2, 'HandleVisibility', 'off'); + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]'); + title('Huddle Test - Geophones') + legend(); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/comp_psd_huddle_test_identification.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:comp_psd_huddle_test_identification +#+caption: Comparison of the PSD of the measured motion during the Huddle test and during the identification +#+RESULTS: +[[file:figs/comp_psd_huddle_test_identification.png]] + +*** Compute transfer functions +#+begin_src matlab + [tf_acc1_est, f] = tfestimate(id.d, id.acc_1, win, [], [], 1/Ts); + [co_acc1_est, ~] = mscohere( id.d, id.acc_1, win, [], [], 1/Ts); + [tf_acc2_est, ~] = tfestimate(id.d, id.acc_2, win, [], [], 1/Ts); + [co_acc2_est, ~] = mscohere( id.d, id.acc_2, win, [], [], 1/Ts); + + [tf_geo1_est, ~] = tfestimate(id.d, id.geo_1, win, [], [], 1/Ts); + [co_geo1_est, ~] = mscohere( id.d, id.geo_1, win, [], [], 1/Ts); + [tf_geo2_est, ~] = tfestimate(id.d, id.geo_2, win, [], [], 1/Ts); + [co_geo2_est, ~] = mscohere( id.d, id.geo_2, win, [], [], 1/Ts); +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, co_acc1_est, '-', 'DisplayName', 'Accelerometer') + set(gca, 'ColorOrderIndex', 1); + plot(f, co_acc2_est, '-', 'HandleVisibility', 'off') + set(gca, 'ColorOrderIndex', 2); + plot(f, co_geo1_est, '-', 'DisplayName', 'Geophone') + set(gca, 'ColorOrderIndex', 2); + plot(f, co_geo2_est, '-', 'HandleVisibility', 'off') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Coherence'); xlabel('Frequency [Hz]'); + hold off; + xlim([2, 2e3]); ylim([0, 1]) + legend(); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/id_sensor_dynamics_coherence.pdf', 'width', 'wide', 'height', 'normal'); +#+end_src + +#+name: fig:id_sensor_dynamics_coherence +#+caption: Coherence for the estimation of the sensor dynamics +#+RESULTS: +[[file:figs/id_sensor_dynamics_coherence.png]] + +Model of the inertial sensors: +#+begin_src matlab + G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)] + G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)] +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_acc1_est./(1i*2*pi*f).^2), '.') + plot(f, abs(tf_acc2_est./(1i*2*pi*f).^2), '.') + plot(f, abs(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $\left[\frac{V}{m/s^2}\right]$'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_acc1_est./(1i*2*pi*f).^2), '.') + plot(f, 180/pi*angle(tf_acc2_est./(1i*2*pi*f).^2), '.') + plot(f, 180/pi*angle(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2], 'x'); + xlim([2, 2e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/id_sensor_dynamics_accelerometers.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:id_sensor_dynamics_accelerometers +#+caption: Identified dynamics of the accelerometers +#+RESULTS: +[[file:figs/id_sensor_dynamics_accelerometers.png]] + + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_geo1_est./(1i*2*pi*f)), '.') + plot(f, abs(tf_geo2_est./(1i*2*pi*f)), '.') + plot(f, abs(squeeze(freqresp(G_geo, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $\left[\frac{V}{m/s}\right]$'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_geo1_est./(1i*2*pi*f)), '.') + plot(f, 180/pi*angle(tf_geo2_est./(1i*2*pi*f)), '.') + plot(f, 180/pi*angle(squeeze(freqresp(G_geo, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2], 'x'); + xlim([0.5, 2e3]); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/id_sensor_dynamics_geophones.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:id_sensor_dynamics_geophones +#+caption: Identified dynamics of the geophones +#+RESULTS: +[[file:figs/id_sensor_dynamics_geophones.png]] + +** Inertial Sensor Noise +*** Load Data +#+begin_src matlab + id = load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab + id.d = detrend(id.d, 0); + id.acc_1 = detrend(id.acc_1, 0); + id.acc_2 = detrend(id.acc_2, 0); + id.geo_1 = detrend(id.geo_1, 0); + id.geo_2 = detrend(id.geo_2, 0); + id.f_meas = detrend(id.f_meas, 0); +#+end_src + +*** ASD of the Measured displacement +#+begin_src matlab + Ts = ht.t(2) - ht.t(1); + win = hann(ceil(10/Ts)); +#+end_src + +#+begin_src matlab + [p_id_d, f] = pwelch(id.d, win, [], [], 1/Ts); + [p_id_acc1, ~] = pwelch(id.acc_1, win, [], [], 1/Ts); + [p_id_acc2, ~] = pwelch(id.acc_2, win, [], [], 1/Ts); + [p_id_geo1, ~] = pwelch(id.geo_1, win, [], [], 1/Ts); + [p_id_geo2, ~] = pwelch(id.geo_2, win, [], [], 1/Ts); + [p_id_fmeas, ~] = pwelch(id.f_meas, win, [], [], 1/Ts); +#+end_src + +Let's use a model of the accelerometer and geophone to compute the motion from the measured voltage. +#+begin_src matlab + G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)] + G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)] +#+end_src + +#+begin_src matlab :exports none + figure; + hold on; + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_id_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'DisplayName', 'Accelerometer'); + set(gca, 'ColorOrderIndex', 1); + plot(f, sqrt(p_id_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_id_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'DisplayName', 'Geophone'); + set(gca, 'ColorOrderIndex', 2); + plot(f, sqrt(p_id_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'HandleVisibility', 'off'); + set(gca, 'ColorOrderIndex', 3); + plot(f, sqrt(p_id_d), 'DisplayName', 'Interferometer'); + hold off; + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]'); + title('Huddle Test') + legend(); +#+end_src + +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/measure_displacement_all_sensors.pdf', 'width', 'full', 'height', 'full'); +#+end_src + +#+name: fig:measure_displacement_all_sensors +#+caption: ASD of the measured displacement as measured by all the sensors +#+RESULTS: +[[file:figs/measure_displacement_all_sensors.png]] + +*** ASD of the Sensor Noise +- [ ] Add formula to estimate the noise from data + +#+begin_src matlab + [coh_acc, ~] = mscohere(id.acc_1, id.acc_2, win, [], [], Fs); + [coh_geo, ~] = mscohere(id.geo_1, id.geo_2, win, [], [], Fs); +#+end_src + +#+begin_src matlab + pN_acc = sqrt(p_id_acc1.*(1 - coh_acc)) .* ... % [V/sqrt(Hz)] + 1./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))); % [m/V] + pN_geo = sqrt(p_id_geo1.*(1 - coh_geo)) .* ... % [V/sqrt(Hz)] + 1./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))); % [m/V] #+end_src #+begin_src matlab :results none figure; hold on; - plot(f, sqrt(p_geo_1)); - plot(f, sqrt(p_geo_2)); + plot(f, sqrt(p_id_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ... + 'DisplayName', 'Accelerometer'); + plot(f, sqrt(p_id_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ... + 'DisplayName', 'Geophone'); + plot(f, pN_acc, '-', 'DisplayName', 'Accelerometers - Noise'); + plot(f, pN_geo, '-', 'DisplayName', 'Geophones - Noise'); hold off; - set(gca, 'xscale', 'log'); - set(gca, 'yscale', 'log'); - xlabel('Frequency [Hz]'); ylabel('ASD Geophones $\left[\frac{m/s}{\sqrt{Hz}}\right]$') + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$'); xlim([1, 5000]); + legend('location', 'northeast'); #+end_src -** Dynamical Uncertainty +#+begin_src matlab :tangle no :exports results :results file replace + exportFig('figs/noise_inertial_sensors_comparison.pdf', 'width', 'full', 'height', 'full'); +#+end_src +#+name: fig:noise_inertial_sensors_comparison +#+caption: Comparison of the computed ASD of the noise of the two inertial sensors +#+RESULTS: +[[file:figs/noise_inertial_sensors_comparison.png]] + +** TODO Inertial Sensor Dynamics Uncertainty +*** Load Data +#+begin_src matlab + id = load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't'); +#+end_src + +#+begin_src matlab + id.d = detrend(id.d, 0); + id.acc_1 = detrend(id.acc_1, 0); + id.acc_2 = detrend(id.acc_2, 0); + id.geo_1 = detrend(id.geo_1, 0); + id.geo_2 = detrend(id.geo_2, 0); + id.f_meas = detrend(id.f_meas, 0); +#+end_src + +*** Compute the dynamics of both sensors +#+begin_src matlab + [tf_acc1_est, f] = tfestimate(id.d, id.acc_1, win, [], [], 1/Ts); + [co_acc1_est, ~] = mscohere( id.d, id.acc_1, win, [], [], 1/Ts); + [tf_acc2_est, ~] = tfestimate(id.d, id.acc_2, win, [], [], 1/Ts); + [co_acc2_est, ~] = mscohere( id.d, id.acc_2, win, [], [], 1/Ts); + + [tf_geo1_est, ~] = tfestimate(id.d, id.geo_1, win, [], [], 1/Ts); + [co_geo1_est, ~] = mscohere( id.d, id.geo_1, win, [], [], 1/Ts); + [tf_geo2_est, ~] = tfestimate(id.d, id.geo_2, win, [], [], 1/Ts); + [co_geo2_est, ~] = mscohere( id.d, id.geo_2, win, [], [], 1/Ts); +#+end_src + +Model of the inertial sensors: +#+begin_src matlab + G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)] + G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)] +#+end_src + +#+begin_src matlab :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(f, abs(tf_acc1_est./(1i*2*pi*f).^2), '.') + plot(f, abs(tf_acc2_est./(1i*2*pi*f).^2), '.') + plot(f, abs(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log'); + ylabel('Amplitude $\left[\frac{V}{m/s^2}\right]$'); set(gca, 'XTickLabel',[]); + hold off; + + ax2 = subplot(2, 1, 2); + hold on; + plot(f, 180/pi*angle(tf_acc1_est./(1i*2*pi*f).^2), '.') + plot(f, 180/pi*angle(tf_acc2_est./(1i*2*pi*f).^2), '.') + plot(f, 180/pi*angle(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-') + set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin'); + ylabel('Phase'); xlabel('Frequency [Hz]'); + hold off; + + linkaxes([ax1,ax2], 'x'); + xlim([2, 2e3]); +#+end_src + +*** Dynamics uncertainty estimation +- [ ] Create a weight for both sensors and a nominal model +- [ ] Make sure the dynamics measured are inside the uncertain model +- [ ] Compare the uncertainties + +*** Dynamical Uncertainty #+begin_src matlab :results none - [T_acc, ~] = tfestimate(acc_1, acc_2, win, [], [], Fs); - [T_geo, ~] = tfestimate(geo_1, geo_2, win, [], [], Fs); + [T_acc, ~] = tfestimate(id.acc_1, id.acc_2, win, [], [], Fs); + [T_geo, ~] = tfestimate(id.geo_1, id.geo_2, win, [], [], Fs); #+end_src #+begin_src matlab :results none :exports none @@ -174,25 +1097,33 @@ Then we compute the Power Spectral Density using =pwelch= function. xlim([1, 5000]); #+end_src -** Sensor Noise +** Compare Time domain Estimation of the displacement +Let's compare the measured accelerations instead of displacement (no integration). + #+begin_src matlab - [coh_acc, ~] = mscohere(acc_1, acc_2, win, [], [], Fs); - [coh_geo, ~] = mscohere(geo_1, geo_2, win, [], [], Fs); + G_lpf = 1/(1 + s/2/pi/5e3); + + acc1_a = lsim(1/G_acc*G_lpf, id.acc_1, id.t); + acc2_a = lsim(1/G_acc*G_lpf, id.acc_2, id.t); #+end_src #+begin_src matlab - pN_acc = p_acc_1.*(1 - coh_acc); - pN_geo = p_geo_1.*(1 - coh_geo); + geo1_a = lsim(1/G_geo*s*G_lpf, id.geo_1, id.t); + geo2_a = lsim(1/G_geo*s*G_lpf, id.geo_2, id.t); #+end_src -#+begin_src matlab :results none +#+begin_src matlab + int_a = lsim(s^2*G_lpf*G_lpf, id.d, id.t); +#+end_src + +#+begin_src matlab figure; hold on; - plot(f, pN_acc, '-', 'DisplayName', 'Accelerometers'); - plot(f, pN_geo, '-', 'DisplayName', 'Geophones'); + plot(id.t, int_a); + plot(id.t, acc1_a); + plot(id.t, acc2_a); + plot(id.t, geo1_a); + plot(id.t, geo2_a); hold off; - set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); - xlabel('Frequency [Hz]'); ylabel('ASD of the Measurement Noise $\left[\frac{m/s}{\sqrt{Hz}}\right]$'); - xlim([1, 5000]); - legend('location', 'northeast'); + xlabel('Time [s]'); ylabel('Acceleration [m]'); #+end_src diff --git a/test-bench.png b/test-bench.png deleted file mode 100644 index b6603f4b26c084e0a74e24a229ab317d26951b21..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 17851 zcmch=zj>0H`xE>)=XJ-@7Xo3*C4I>J67medFFE~mpnfKvkNpB+K2DIpz(6S% z4_99}%*#p2)5ke?UF|Xi!VkHteakc`Z*4j-k9D{8Y-@bN4$212{uszcYi3I63VFyv z_u3)z$ql3pgDy)?_kmnMiil$NmCGse`V~;2B0D-He<^L^Kao~RgGsw#<#tUCx#s0O zb!&$Kw_r^^M&7RBu#v&Q3j@I;4TIc#OxoJoEq`cqA>yr7@i$tPqB8DlyvfYGRE-xD z9DUq+vNyTpebpnEV0 zLBV;XQnfh0`^Y=>(1Z2n4P31=qp+~>6%mnFuU_5O)uo&i`+9EGbRnC_Z5AmwZTg!1~nJ zmf#%@)c^?!YWgSt782rrYQuX+uV|bd`)=>-a4LG(K$+tD`}-ZI8oi54N)$uF^$iTp zQO9)m{0%YHZYjAh=;GJZ)I`7`V^VokM?V_mTplnVSGTpY()aL?LL!mipp4hAIqKD! zXc3WT5%UWRcMS|E9A?|XFR`;{ynDyH_aOJxtBdF9=*X^$iwkLL7QBA_I`KLc6;)~Hh{y}0e3gd}9^91SxN_}Unf^=SNkn2|VtiU!wXr61`rE6?AHTa# zqgdZ#`BO!}#Tn{&Q@(31<=oJ;ii`0-=Q&amnuCnwp*gB)y9 zHhs|)G_;iD-QC?@Q{I^GuuOdYrVJ&Jg;xK826w$ zVrXhTrd1T&aA*)XTZdt+)sv0hLLdl?nUeYhe!NEBKy7dPK%RBHE!b-a-g(R>Y1IuL zNE|B+2AWFtW*V^B+@-YdBu!x?8}86 zc-7R@^!4>Y-iH(S{P5V6XV1t%XjyK(Eh-WW*q9mmR;Yy!<7Q^g>4aXkH8#F@vq#C@ ze{W-Ug@yd_-n4%C#FAQ@P(%-)Qb*14^!X5kvnQYzVRMrg~@Qznnvt2QZD(f(9rw7R_?+e9MKP@VOv@Y6S`yHDu^Gd69QC^^ z&SL855@%cVP5!5NMJ8I)W+|(v=E!htb zRh6~Wm$2ftaN?`XOw!z0HU~(?P&16G4+#3;E`dcHvb)-YL&C3Jy*k&rKLcML&WFLQ zjEy_1Y@6@s>z@x*@VBwy28+RMwk`iAKN^g^WxX4x!#;B}Y?a#ii%e?$kbCU%;vy40 zHGr2Y+@I|7^&&F2P8>R#LfLB$GFk5~q^e5%Hp9f6&V<1MG^epqDP)iE_4zrLgOv|y(AE%n z;@;ZAV#a-dWnjE7DU|YMSN$IQd4q zxNs#ydyz@_$`u)x(JR&E0)mzCR=t9Pg1rB8evpkt#>dBF<{=uV2MeWX%x>p7IXU0) zG`A%tCEdMqr*x~Fs9n5-YYlC1#z={a8$f;qo0Dysa^D0)7GpP6g!EqMb@1?5w<#F0 zIN4v7hezL>1OMMo7+IcN_qWxydV7=K%kkd5L5rKwvbnoTe};6 z9?s0noOUhw|Lq^3dE=~1%9IpQ_Z;{T7dAt!d6?5*6lQx5jPhJ1dxl0_oz(QMS2ym` zcoltnPE!D{L1j*Mm%(I61y4R$t<@xp&OzhRe9^*a4+q@^v&xAbGJ7(E?>cpNl1gMe zEbA$1^V`p^^zVu0FS*DQ7HFXSv+Olh>fcg7orNFBZt6=q9YOdYwNcDyZHKFScH9%C zb*?jEH{b5)ou%nmDV>OV1fQ--j3m^!E4)mj!hL&c7pZ|QI`Zs82}R$6CN*6%q_C^{ zIB0g0X0GX_e{SCRslCuxxiGB*yQo&QaT?$yyTvjd!EKam4t$+Lg7r%#6F#X}m38vG zTS4U<@9D2TMXO+SEcnTAHU>=*!!_aEM`P$D2HP^WPn1{pP!|P_BQLu?eZx!E3$I-8 z6sm899FxB$pW;KV2|BQ6MX#ag;;B43H0H=YkuOB)ArWoVIDMy=#5{OjV(jYeIhWNL z{c;0#dGz%;J@hQ1>hP`A#<=P(YP?Jiw{hSqwx>vAR;Rq*Ar&x+@Lm&*diHY7g1i;6uYB%+Q|2=6oEjw|9yj}NSlL-= za;Bpf7lUgbw;i7K?6&H!GOeIlyjn7KeU6h~vG|^JiH1F_oly_gexD}~K6pO3TDv&C z%xUU+>+^7gU3@t1#pPv*G>lb&6yK&Mv@}bn<&T&|-|hUB3NDJUcq$CU*;-J6j9}+H z$vtO2XSl57|1d}DX%tOV9r`a?f6hvRP;^M4&>`-ZUB}nRLWln570gYCRE$bP(!IuL z)iYv-&*E@Rf0QxpEUL=~FKc9Rv3)g*Be~ygfIFojy6WuKYzGS@iRSA4mO9w+s@D+51)vc*t+Y|6sNx{auaMh^5ch|(9jJ8B40+{ za4cwWtbDTH+RTU#R=Rl5X&J~M!%i1~FyGR7&M0~bvLv8ML{t&UD-=6gSB`LHXSKL* zy8Ojbnw=pJj+B=&kZ)FhAn6pDUC5tuHrS*&@O`#;UO|*)W3zA&xFPvyatzw=H~K2E zpXP$DC5;Rtog)Ha&`BOCuRuY(vz#`$s`l}Mf3+)novZYwiNqbd?#_R*a}5ty^VpYm zgPHNu**mN^BbPTq5!k|kR*vAzetv3veQ#MB{;34|d|8<06y5Z<{jTZdU{h*YXgTY9 zJ#C)*7^^ymGwl7?0`siM1~~Fvvf{aoldD~w$DvTV$#6+r>|~0-QW9ckEg}8}XI--O zi*;P3S|JwBS6o6C9-RUmpFcrU6J~yAAh2F~bC2e7#t7|Eo9m{5JkWBnpAP&zhb$9y z>bd>0<8SYES=q%sRc3qSNpVK*f$>=O!Q7@xI90X16Etw}e&<&%GA^!hVfy}qiBXAo zC~6AAy?8uDMOT%t9jSR6^KUv_$L^r<7RN8ca5?PIv3*JE2wZaIgqcXkFfEj~#P-=o zw4Z@WqqIo9dG>T_fFLs{cx84{G}1J`VwV z6TC-H=QcZ(eMbm=MQr36(~^AQlvm(v+V0QOh@bT%-yHH|QqJrkHsn*!ZR5G`I*Dug zM1CEIku$$ZAG7bt~T)S8Js`QvYP*_jG*5f^6)e`q-lko)WstvQ+86MBnyfN`OHfWD07V78n1!Wv>z&=KR=u-0>eqc<=04r!g=baiv{65At-P|K3@b=_9lgB{*b z_VN5MQ(A?nXG!d{`g*HrDoakBRbKXIVBC{D(~Y6_&eVSD7mG zjL}i9wxIZPER^=_a&vs+%g+rSm4T`I1`x#%rIip!Y=q@X>swQ=_;1oEmgr`|8oS?x`S30 z+VFcCIv(l9!k((hLdi@6+x4B#}5O`n7xRki(GJ9{2T2p%|*E!BFtybSIBQUSFrBTS8B$0Sbdfh zi0?r87uL917v-i+8}BkHT-l^o_-d$qA-c2iZ(eeWfJVG}ivU;X`{(P^Zttd$KR<0; ze;VA6B)*&$NE*Yf1^9>G_W9yYBJY|EbdB915)-q8I)FW*MhbGx8-xgqPKW|9jo_>NXscY#5uct zrjv-r%--PJMtp?&w=gK(0fn1e5F{Yi@Qd?pgPmPg?Cw$3IIBzHdBiLcjsLo{*JL&-9(u?zKFbNs`J;!UFHUAc=Q7aU3v@k;pb&87cVJY4$!sw)1E{xD=+0Z4OtZ;1?>%R3|T zHu85{D9`{qxYTQGuLvP?Fk*0u_?vn7C4F_%_dTsrZTj9u>^%qi2Lyi|?cy!a!)n`x z6x3D>-zvEmAv7kyaEH>E%jGOijJ1*W0H(h8yUjwrwOVjgA9QZoRxe%-_P=Vv7fGW+ zE`4M4lyU4kx;^UFE>17H)~$k$gV{xwbXos%R)N-V>Vat8sr^EIzib`wg~D2iGG|-u zVygHZ6#3{ptBwR?-Mx#OCKd8;AQfm4$))F+E;PqvP)0TCV|pdt13Cno5}guT5Zys{ z@X*Jbx|s{YkZSUQ`~rA1AB^RbaTEi(JvxQCx`_uzu0}}@k%EMf4@5;moQx*EqWzXP zG7xM;vbpq^C%cqVWFcg#Q5I3vLaUs~#(~T(AAerYe(~x#)%`D;E=F(qqW-4DehY@) z)S%<0s0hfeQkn~qcpjPJZ(y(Of_&d$p|L#wIuByAjW>ur z===(GMK}2jt5?4by$y93b+XsqVU!8qxD3DYVwg$ZcymOuEv=td>I5T{GA}Lui{(>J z1Fa@XrRKNp31bzJY_4zK=`2!jo>b~-l0b(3zNmJT;`o49iMCu*& zI^CgiR9$8*l5<&blp2)ecAN_m4cZg+6P6dK{wPOELlnqDi!uYd8fE(HJwV|Rn#V_L zsd?TtYuxWxfEm*#+A_fT{-@;&I~=O;vLzd#Q;>sZ8iPMz#*z>DP(;my;W=A03yRg&d`Fe==%}PZ!hM2%E^)GE?MGZF7{!0S*`Ftl% zvlOzn9@IxoG9zk@_#k!Dcg?41iC@UTG=$eBL3K0OJ`RRBPe|(@+70^{_Msk0n0~l> z{u;|FxtrS2I8Mp0ltG@m-B}(}?V+HAYitZ<{!Ja`l3C-6;9Mwq5bXU@du-0gJ9}v9 zmWbeJyx*AQZWbG|qU_!ij*%`@e2Jr6l`Rsj49)o$neQ(Py=oSi{KS8S`j|DBY0n{Z zyn<2#RDF%RfZ;bg6WOb3?IKaQ0>u2CmA%1Q}0wJ~?gPma?aV^k(zRE}A4{acxumz8*2 z-9Nz_WWtlCd3OaAU~be3bVvGgN^V8jsv11_7KIw|CKSBdH=7-1i$8~}SlvKOOyp2Z zB5t0(sZFMeU~Hj~_*A1)0RV=Aw7kos@GIt+_|Jq%STy{*L$}#WRyn_XOnY z{>FW>3tq}90Nu$bWJqf5Q+2AB67PCe`_DElBrDK%|e zGE~Hnp03{!p~iMR{WIsLAbCqaooNwVjWL^^fft!wOu==8r5*o%975MU#Ck$f1wtN- z>qKB9JkJXW+4NjMF% zmFIQWRVhj+?BSIoNNDxF>)7)~U2R=;=tJVxO-PWUV*yXJNtEfF>YPxCO=(|2d2N@a z9l6>C0^@|a`H*WWB-O7z65_@%kxh9F=_Awdk9bG#*U?-2!PbGEycZ51$48G8{cLCS zB7ga8S7eYvU$FleLJSOIAqM6p4>3U3+(sZGWNXoDI3nAjD!(AFv7-OvMZfHHzJoKy zpmPVjY<+y>4UV^t$k^DV_!O{5@-mk==)T>!M%N72L=V}kSm*lo62lA;uge9BD4^QmdsU5f(dee*^r-3i4#h5kM7zJ zve0%9|9Vfz>3SZ99l8>SqH4yioac^f9T|ypOCvjFr`v)V@kFGx(&*_D@ZNK8Uit=; z6C23z>-N#McmDa-hsT|6B;62y%1LJ$V6VwA#h8~;cgjm9g|YH1P7=?o>uc9JW_ z$9{7Yv+i5h*yw(&$nGfFtPZ7RBio?ROe#&9bX4U)M!fYUSPqX!tN0(Wnn_cx3dU1~ zwbh2v=QjE)vK4i_tVx4Q4m>oie=$kHNVcIw^co4|@M7`0jJQPi*WYNmz^mU6GgDm~ zJ)ZABLPl{rH8_9gExgpt?2MIePZ5#Czu(%A3AyCa6LP^ZdTF?xh@?y4$}2cqCp1*4 z(+ZAK`uR&dI4q3dVdL;nT*C~HG9(1e={m?}-J;(-x7Y0Fkr{rkwT-EpaTg=v9Y|$AVj+BM*slwz|RFyP0 zx1(3w&%zm4BrUleH|1}G_bOwtg)1)izJ`wzU>`_c&8EV!J?m0TnBrH=h+ zQY|;NO-x-f3zF?=czA4`@%*5?e&`Fs@0 z+_Ca#0J(q6EAh{^3#3RXecNg8IavYx`6wT^FwIdGqX`Cd29zBYe?tx|+8*~eI;YJR zim+=w&0K1}2y8LYSqU4`Ca@^KH0$tePFa)vHS zZEO~dT{@km869$cxHCjF;zj;~#D$F$=5*a684v3z#?7DcyB_6|0*h z-_Wr&7O^fRD_rA}mi4O>xlM5po(VU-U8wImt!|n;Vv2}JLdvdODw(?MM*S~H ztZUp~4BHu@{>WB!xem$flyUA%aU#|5ov*!+;Lvjbd1X?`(*;7B?}Y}4NJvtbMJq;(Bz`!FE-ncyRYda^y zwohB~DgO!nX32MSGuavGJvR%R?s39Vs7eL)qZHtkOHbItdeHX>Y@lYs{xWv!geKP( z7z!-e_xAdfz$N=)ag?9s^21wbk1s)E{r4ku#zb9pP%QWh8f3Y3dGYE&3%U~TPXwHt z$z745v3%Y$If2Lr0nO_yc?-I^vZwjJ!|P`aPPiVl=v~rz5eV2zNCB?;a4TPS1^D}XT|Bhwv4P2jU+Uwl)!|((~g*U~5 zT@o>Hq@BpU*J!XYcKJl9Abii##m;UaoMGpo9*>TRjp;#Hn}ZC;+O_wBquDNl|IfBn zyzMObZ)?t7O~WG~xK(c+?PBjM6nDqmTs;cnMqk`$uy0E^ckUcx;M%Xc?xdKF$p-1d z_wVP@Z}M-tN~cN-rmTm|Wv?n6Ip$946@6hHT33`-Rc&7Epy6{{?m_ZJ9KarrysNAV zO*faQoL#Ar!n|YQYijc%UJeM~ZecgO=n?^h6r-=hklx)QS1c^~SPE8dxm#tK&e1zC zZSrn5HH^7fk_z3SA1*if9E3Kj8qCVZ?xI<4^2$(_*LEH-l^U397Y84aZN}GDRIiFu z%0It1cg^9t7o`{ZrmhfgwepkAeG}3()Rt5ST0#gy%m>n&L9MPPt!Fe6JR-Au+LYOora zdAIy=>twzabO4_bjhg9sMSgHQV%?%^%a*5<$o(<%@}BqATEeSx<;S*GU)D<7#8u8+ zrRjef6i)vw9p&UE@v>9~>@@yPN|;zmN)cCjuDJ$=!CrfsjlKOBe#eFSc(_K*@g!#$ z1+5tHqr-;G(B$zMuLq(nhG!$GD=;Vvt&&%KRAKZLkT$Be4y%Yl4gb|w^v85pAs zTTX#%cZW{{$*1$jd&m_X86Uo-A&~P5t342f^q^1B>BCI~%;U!~vGcQ0hPW&B{%Ae) zcV_?NBohov^a3(Yk(Hb#%3j64odvg6F}Ylk_K(G(Zs%wSt$LHu*_yE=LR59Re&+D7 z=H++^B+`HQP-{7;fhA(gNA|&an}i{aS&EA9ZcmwAuJWMX1x4Hrv+6Bu`J|8kK$znV zh&XmDcu{U%zR=a_kX_qCNmqo4igd^**^Z8jYKZ=bu0``jL9(XH7ZB(m+5-)2VkqkY zqYQ*jaz98G3SlQd!xIAKH&xg&4Wi174#S7Hug&>Ml;)6qKyqS3qi0o}^pRfcD};S! ziErnlf|ixmkekzMc^hz?(t+%7gsm*}TMA@7V0G&f+ds`tDcWIj6;;fk(B_pa+XK0% zc8X`3h}v6hbGTuVGv8}ZAXhGoHyGEc>lw#mGO)@Z6`RgaCpVRwERhUS-7 zjX9F;#M~u(%X)sna`#u9v`BsQ(;rV+l?KinB@M@SaF0H_IUr4j8Gbs7X^rk8YS%*~ ziYcHn1?P^rU4O^&R#UJ+3(rUMW1YBAXLq5l#|@^`f5xHW#}?=zC{y#T04X+Xp}HVD ztT2$8t=r*N06SK!os~}R_8!YdOK>yeVO3+P`o z4rcbff*nAND39}-&5k4?vkxbrOfj}NH_FGxoV^~lk>}}qL&&W#;l0vr2ZBJzViWW` z_R5tjLAE^1WhKz~X_#acwZCey=JYpJ_mK6tLpU?uH{Yp5h46uoEg~1-eUt5=P!ab= zLEHmOZT?aMPV1)41eDz`(!s2+ZJC+38)A~Tq{xgPhgbg0(84yc)8A*(K=v~X=Ea9M z;#E}!ejwO_?{|*)Cf9P`a>|MHpYnaDd$9I|2PwO&UBSAJ-POtCe_?*Hm!!0?k35mj zu?s01H)BqxWO{?q)jBsCi+w|+PZd|G?2J9fQsXzvzCzvNs&%4F@wCws< zI`o~|W5{Raxm>4FMCpL+LX7WjKz3-6+av{jSMRf=+k3vHQ%IXYn*(?!A7Z{uF-Z&C z^Kkd%dB5*zVE5hxE_9_W?cqq1+4HU7<$TPzedG|(T=eySwr776A0$oNaWA*@$2XB1 zD?^1T;K#g`b=!D_VV7E zoy0+ro3C}FrRkdRU3BLVL;dP1$C2IGmx$`xsk>_%dTQUDjhmNSUiW+KKAD2Y<@=V? zE-WrSszR-bd9?1D_onfU^rZ{kH8K(bQmdjOL7*=9i>>79dcokjqSd^=$upXqg%xC? z{ZgPF8gLlL7s2s_0AE^Myh656c2w|LfU>%E@7HK!W23{G4nFw3gA6nNViG~YmBH~s zh2%QrZqAs!9wRp`TUWaBP-d zlrO|2wdnLS%=f8Od1q&oUhe(A$Fo-d`4D%dQucYAx1L+Y{(=fYwDRq4`qH%2_FZ)A z!_|vc%R!&p3f47@y!Mz6Q-tc?k)G1B_OYz4XAe3D~$+v~%FFdD&z!!X(0XRLvHUo7ily)S?7Qu@*| z`HYJLnI=e3_&9(}&}r}OeW(h{__B?>dq16+(h^@GMR*L$_!^N^cWGG#pGLh78Se1F ze7=^4qg3kcS_1fA1XLF*Fnw5Qr|m(rjH+ zQIQP9P21bs95N0!WB3#m6tI!|(}AKXT#8j4`;rg;y82tCgTu12GB1sw*FCd7Kc#!M zY873$r04)UU~&7i&PGt+)7<1DU~A}*Z(3=Pvbc8g{hz;s3NeXO^*qg$e>oj4D7bu$ zOX0J>=Y%ld56MHlTP47+Wi{CFFx(RjPUKZE5p&)_2-5k{3OD;jpxGu`XA*5eIBk_c7doW#Vhh{$ns=HAgp^9pFk`gAjpSCN+EP5zrer1cf3 zErIx~yVsnj_hddvdS$j zNH-w4oDJJyh-fe9?7R(VEugC(0yk!oFUuz$9Bgjnwr{i%gXR-NUn~25>V&p1Nglu7 zrsG_-YGXQRw@NjccI{Zts5Fw@Gp$%XK;EnzzV$S6F8~r++V|>mKcV7r_#%~oJHDe+ z?kLFlW3?;J{8>t%g=kLyS36d#p#>R|W^ugzv}er1eyII)w>ojxv;9P&AhN0&M6u#b zP>gt0 zvBEObHG16Mlz~&;dOb(3RZTBketeSk{c{gIpUD<+VZB!B7T*v&O}&?fFAPm5-QoqU zW>2sGVD=^1K6Fi9C2)PvmOaK#8e;h7@NXfq|-Hzm%>}4VF=F3RlzW=kk^DYWtrUSis7+Xa=q#^*No@qEFd{ z&6Nm5z2}qice5ML`|EGl)2r1aP4 z6+@N%7#+m{y`|O*`UF%sr~)26`gv-TT4~*Ik)1sdHL{Gib%nieX;Ewaa%0=>5Gz_E z%=vWatvo1a5TT9jhkv*{n%5ML_BTfCPq%XtcMlH_)s9v^3>Ih>0I{kg9DYqgqPVm) zJ~^4*E_~MvqH(y$vfh3ck*YLL>w)_F39%q)i^|0Ch9-}dukoHloE*vb1p`GVYgD?B z9u`#)OUJ(bO^XrvctqsrSN1sEJ|le@jbCB!*S|@#7Yuc7&w!q+nD-FQjP>*N$j#u{{Ig%?c>; z8G%4Rz}tW<-6ry6_jo<9eIc5L5KI|3?QhXK;sj>>KjhBN&c0Yn=P~FWICddXylO60 zb$u>bZrEX|Kg%ZM??;k)dAXn=ZvD~Y$3FPV2A>{Q3!r%%9(*H5#&M%>9db=fY`!Q- zs;H$UucRR2Eub_$4=tBaPH|f2ypN17PYXIQ>^rmrX}WPMYVsY*#}PCiG{CfBpIuHc{i`xB4?7 z*R$;n7hNAjQS|A* zR^rt`#co=0w&oc*;5(uocxM?2RMDrW`|TMoU!s9v{k3YK6p$YOxgU9Cl2YjmG6#(h zajiklDKWfo5Gbb)K>+lDBQZcC`S)iT3y@?!di02ikMF*Zk586p`^iqh+2W*U`x3ys zkJb$y5S#Gbv18)w*&h%b%l9LnW-+c@OuzgkQHB~|t_bN;yW{rG0I?UBmw$XTU{6*4 z!y(mnFrSFjQ4af{n==3J$9?CHXHjU6Q@D{Fxc+4lGR;rTWw&`4 zmYMl{)@0*_*AtG&%?_GNY;5L#b@vW0oLtL(^X3ix zcFDtjpq}wvV%SUS)A<%CA>po|Vb$1kHukWvWY4H`whWvXnKR*xs@`M>kfcMF`tS?Z>;JaQHlIf|U=v(P#&br^*u+HK>PI>|b%z6@f5QeW z6d*zdGxDw6q8i~2sEj?1D)~m@@&Fy$7h=r3S(K z$7H0W4HER%3Fm^LANMs4fhU1clSv;yR&fqhg&0;GiYaC-QVT^0S+D?PCV8@N6qte9abTc;86eEJd)Vye*4EUwYqVi&AAS%0&n_nO}@^+O)?`U zCZ56vHcD^V zOyebk89Nm)>D~?coI68`a^CaW)#Ve8X6%jtq-=b?3Qzfu+u8uEZ2_?Y+%6WA zuZ97G`}M^IaGSW}<1k=Sxnp1e|MvbC>D#Og*L76w+}tE0iIcO9|3(eofL%ne@Es4g zm+lp4evK8E527+D+23I&pFXNGfM$s+Mhiv=l#B%mczXHn)LoKyk{8qBUf9zC`on8keHSV$<@wg>IJn#rEwHBr{BM zaIWEV3G6O_cYW}jZ3CmutPC^HsE<7=SM85?G!jTGF(OwN?h8oB(1PIB~oebV;y|;fox! z_H5NDG0NwVdrW^{xTWkd-Rb9U3?^Zc?g#5)oa&aU_z2P<^CWZ$0gkv7=PNy!4EtHnrDY2Hoimh z;4d#G0UJT^ZA)iS({nsntJDNiC$q;|fh*q%`%(Kdp=mU{K~c`{faF^%Pe&U(?T_}>xax zF<#*Fyi)V`bC%C78Ug}^^AU?rfUm0(Y$KCdK{?pj*>`|v?A7bncR*~qxVStVVBl4I z4S2vi5WWwLjc?z&MFu9osXLz8&<5_U4w5N7mx9D9{O)v|;wLjF2c*Ru#9ccca`9+SePQk!)sO??WMsff6^`N5I&;){NAsbsnkO7Ty4x7(+VcuhH@#?HVaB zU788sx20y}y$39ZqylFqbVC?qaFFur#_jO_@}x_HyQDX5iaP*QL=-HEEOCqLek%{u zFPLa+zpXT9>HQA)#T^0QX&kaE=a8#13=UQX7CSo&Ep6>9pUuBtzX0eQkx(1<@4FtT z&@a+cgM}mTU%O2J`$;r6$3CKI+>Di(xvU5z=JM*ce!zimfye{3Z|Cw3*`Cv`j*bp+ z37J-~Ff+S&Sm|&`+qdO`%E;^NWFN2U4FZ* zj#*HWA6db`fx+OwG!V$ZuC@&#Jj51UGF=N8EBryg9*k8-fQ=h+AmIR(wu>zAaZ>L< z+4m30I(U(Vh2%&40ICfD;iSD7faa|gZSnDUOeMvVGcq!^h$r$9f4iPh(G^u`u9kkd zTy6!%dUv#_lY_*pA%yymz<7zkdY3xLK)nk96v={dFtjxTx1i$U4Vwa#MI~zoSo3Nt zz^K2i1@~_NW>whLTM`}xa-WH*Y4=M3trP;4S}95r6Y&QE zz#EE^Qppt2kcKHAI+_dIb0A?Z0#aX!s|OWLRol^1;MWw#g383#4Tp$s-e|K<>%PRC z{a01k76Ratp5%qn9*YDe) zd$<7``j7DRhWCQx3d}+(67giOQn<5cmx3hI0z;8+v_EhhjUsYz0((qGo^NQPSM ze5VrFp9HM?MI{||Obf_;5K=}2iS>p12552sXfZU`BhAqF`A*=}?2xc;QfCo7$5($swY+!4Y2P%NR!~fxIxnR2r%*ZXuSOs_z7FWR@;xHD9-c& z7bPeX^YrRSk!qPl504^M_U^xK02^8P3GV`6M5Jn=*5!vtO>J#%;|mY~{y!YMaAvkl zF@}5n!^7rGz*<;~o&OK5<`;$keeY8naKk!;g@vJ$Vra;tISNJ^sIF3P5>M9PaCo!d z-%;Teo$w3XN^{~85^79bT)7K9TRPp=9S&X(?%hikErG1Z!4BL^BB->u zgBR0Dk!Uwom_uZ5pB(||77}J>S`79kzX0sg0(+S>r0u$)XQJva=-4E#N=i#hZ_eZ9 zff0Ie@E*Y8)+#!PddR;`)OkTcLn&YlytG+tCN%?$aU=)Z|B=@8e{F{QzjP{zVf^_! z_)GNE6YY_jv9Yn;Jw55bg$+t~S=ph5BWVDLw3P8lKO)JPx87+a6yE|kQv;S0u+x&< zc3|)d0WN?F3qta}4iA31pz)r9- zpfWC?4(_Q3Nv938J4l}0r%yjG!}4Cfya{F*YCy6{gAz4NgAWj+ZA5_JC>Cr6-vUrb zDt|6sy!cuMM!vPR1)BtA8YjpLV8azAu$BdiZK4^Lm5r_D$#=mQ^qiyrSTCk47{VkN z9KXI22F_QR00K_z7WY%o%>ZS%JZxxHiis8iLI&?`2%PEmq8&`ydzqf`40?xuyp9V@y z(lY@2sSa=qPXOOBNTwunW-G2Jk zL8cGb2_J5~1w9`yKp4;;B?H8h4h)EKP|}@N`Lwdi(XXbHf)EIm4C%c9MPPdZ)S!RA z6;fyvY`nB@2aKVOxDM7muA@jN3D+Fx*#li4&Q}cvT@6q>)AMQM0hEx7JPiYncDX(c z*p*3BYh-kRFKm+xKVB`Wto#PRKFI_Q?x3sR6F7bS@mRBC!afnqK3VCC{sw@fp1t{b zM^Zwqcbn)0U0Dzck3r)G^sCNZj#y}4c?I;8o*j8tOr@p|=VU7W9Eo0E>1Fa=-tG!v`%s#XwNEU!I=;QEWo08&z7xi+~WTK>` zWF>M648!4O7vo)h{dqmjqlY~qicB;Nmra*676%6h#clz@st@uxxM|i(%RoNRu(&`s z4R%SwH0*n*v{qyR?BM-oeL*;o=2T8D1K8iqD*!;ykQzJy-^$oykw2}itS&=9|BN!G z<7xBf>1^M{*GR5n&^&o)zYMzcD1UL1cUV!e!I2y)!%q`zi}NyTuxXrr-dd zk{1NAIo%TMi>n9d^upg6n0hae&I{NVC@MfWBLG@G-4F=qsz$97p4>u#R|}998E^`7 zTz^Ch-1<&`M?ZogU47YG@}VwDtBYIA*+g%B{6rbt^M_QG&%>BY%w^(+bFV^@5)xp* zPaiduoJy~p!_}D%1GC(fU8r&L&jhT6*Gk@k0POey%oW)@L8I_fRd5UN4gz5Bo@kie zsID1iSm>Jr{bnkls1p=qAg#Ghp7NogAt$i%9#Tv&G=VM{V%H-|m)gbBl9Ps!_X?Pa zBLB;@l%1#)*U31d?_%o^aCIM?oF+nE^k@4#Iswd?R5OCGApEKmcbf$5r-`=k&@5k8 zX5~w(pR%dJ#*Py32_GKetHHh(11R5$K@eFUa@Qh(Bs?>o7=-{12Lo`LkLl9DN4`%J-O z>JI3Lf;|P`Z(jhH-NV7(D0=de%0#_1UU260oE$w+l>96SPKf7DgwA1sA>Q#&WRqKt1^qgeYLIpgY7`%?RGwfWcT0H2uPJq-I@O4+)M}sRr!CMFbOER!0 z!v@*1yM$=)D>?ONCVI;P1{vbncKA5L&Ti&@v|?*j(hc<#&j0aBkN+#Ki~r3`I1f30 z{yZ2VyK%sYhqj3)Pr2)08h|HAbI=ED0>6f!yA6{P+3O8|o#MYIk;?5lQmf*kVc^9p Okh?lY+Euq5p8Xey9UVLX diff --git a/test-bench.svg b/test-bench.svg deleted file mode 100644 index 23b1810..0000000 --- a/test-bench.svg +++ /dev/null @@ -1,312 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - Rigid Plate - Shaker - HorizontalInertial Sensors - Granite - Wire / Suspension - - - Gantry - - - Interferometer - -