Update all figures

2020-09-22 09:51:26 +02:00
parent e63f256600
commit 5303413aa0
26 changed files with 1728 additions and 1567 deletions
--- a/matlab/index.html
+++ b/matlab/index.html
--- a/matlab/index.org
+++ b/matlab/index.org
@@ -75,7 +75,7 @@ Let's consider the sensor fusion architecture shown on figure [[fig:fusion_two_n
 $\tilde{n}_1$ and $\tilde{n}_2$ are normalized white noise:
 #+name: eq:normalized_noise
 \begin{equation}
-  \Phi_{\tilde{n}_1}(\omega) = \Phi_{\tilde{n}_1}(\omega) = 1
+  \Phi_{\tilde{n}_1}(\omega) = \Phi_{\tilde{n}_2}(\omega) = 1
 \end{equation}

 #+name: fig:fusion_two_noisy_sensors_weights
@@ -1027,7 +1027,7 @@ If $H_1(s)$ and $H_2(s)$ are designed such that
 The maximum phase shift due to dynamic uncertainty at frequency $\omega$ will be $\Delta\phi_\text{max}(\omega)$.

 ** Requirements as an $\mathcal{H}_\infty$ norm
-We know try to express this requirement in terms of an $\mathcal{H}_\infty$ norm.
+We now try to express this requirement in terms of an $\mathcal{H}_\infty$ norm.

 Let's define one weight $w_\phi(s)$ that represents the maximum wanted phase uncertainty:
 \[ |w_{\phi}(j\omega)|^{-1} \approx \sin(\Delta\phi_{\text{max}}(\omega)), \quad \forall\omega \]
@@ -1738,7 +1738,7 @@ We define the generalized plant that will be used for the mixed synthesis.
 #+end_src

 ** Mixed $\mathcal{H}_2$ / $\mathcal{H}_\infty$ Synthesis
-The mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis is performed below.
+The mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis is performed below.
 #+begin_src matlab
  Nmeas = 1; Ncon = 1; Nz2 = 2;

@@ -1793,7 +1793,7 @@ The obtained complementary filters are shown in Fig. [[fig:comp_filters_mixed_sy
 #+end_src

 #+NAME: fig:comp_filters_mixed_synthesis
-#+CAPTION: Obtained complementary filters after mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/comp_filters_mixed_synthesis.png][png]], [[./figs/comp_filters_mixed_synthesis.pdf][pdf]])
+#+CAPTION: Obtained complementary filters after mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis ([[./figs/comp_filters_mixed_synthesis.png][png]], [[./figs/comp_filters_mixed_synthesis.pdf][pdf]])
 [[file:figs/comp_filters_mixed_synthesis.png]]

 ** Obtained Super Sensor's noise
@@ -1824,7 +1824,7 @@ The PSD and CPS of the super sensor's noise are shown in Fig. [[fig:psd_super_se
 #+end_src

 #+NAME: fig:psd_super_sensor_mixed_syn
-#+CAPTION: Power Spectral Density of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/psd_super_sensor_mixed_syn.png][png]], [[./figs/psd_super_sensor_mixed_syn.pdf][pdf]])
+#+CAPTION: Power Spectral Density of the Super Sensor obtained with the mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis ([[./figs/psd_super_sensor_mixed_syn.png][png]], [[./figs/psd_super_sensor_mixed_syn.pdf][pdf]])
 [[file:figs/psd_super_sensor_mixed_syn.png]]


@@ -1854,7 +1854,7 @@ The PSD and CPS of the super sensor's noise are shown in Fig. [[fig:psd_super_se
 #+end_src

 #+NAME: fig:cps_super_sensor_mixed_syn
-#+CAPTION: Cumulative Power Spectrum of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/cps_super_sensor_mixed_syn.png][png]], [[./figs/cps_super_sensor_mixed_syn.pdf][pdf]])
+#+CAPTION: Cumulative Power Spectrum of the Super Sensor obtained with the mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis ([[./figs/cps_super_sensor_mixed_syn.png][png]], [[./figs/cps_super_sensor_mixed_syn.pdf][pdf]])
 [[file:figs/cps_super_sensor_mixed_syn.png]]

 ** Obtained Super Sensor's Uncertainty
@@ -2285,7 +2285,7 @@ The PSD and CPS of the super sensor's noise obtained with the CVX toolbox and =h
 #+end_src

 #+NAME: fig:psd_compare_cvx_h2i
-#+CAPTION: Power Spectral Density of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/psd_compare_cvx_h2i.png][png]], [[./figs/psd_compare_cvx_h2i.pdf][pdf]])
+#+CAPTION: Power Spectral Density of the Super Sensor obtained with the mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis ([[./figs/psd_compare_cvx_h2i.png][png]], [[./figs/psd_compare_cvx_h2i.pdf][pdf]])
 [[file:figs/psd_compare_cvx_h2i.png]]


@@ -2317,7 +2317,7 @@ The PSD and CPS of the super sensor's noise obtained with the CVX toolbox and =h
 #+end_src

 #+NAME: fig:cps_compare_cvx_h2i
-#+CAPTION: Cumulative Power Spectrum of the Super Sensor obtained with the mixed $\mathcal{H}_2$/$\mathcal{H}_\infty$ synthesis ([[./figs/cps_compare_cvx_h2i.png][png]], [[./figs/cps_compare_cvx_h2i.pdf][pdf]])
+#+CAPTION: Cumulative Power Spectrum of the Super Sensor obtained with the mixed $\mathcal{H}_2/\mathcal{H}_\infty$ synthesis ([[./figs/cps_compare_cvx_h2i.png][png]], [[./figs/cps_compare_cvx_h2i.pdf][pdf]])
 [[file:figs/cps_compare_cvx_h2i.png]]

 ** Obtained Super Sensor's Uncertainty
@@ -3509,7 +3509,7 @@ Since the input signal $x$ and the instrumental noise $n$ are incoherent:
  |\Phi_{\hat{x}}(\omega)| = |\Phi_n(\omega)| + |\Phi_x(\omega)|
 \end{equation}

-From equations [[eq:coh_bis]] and [[eq:incoherent_noise]], we finally obtain
+From equations eqref:eq:coh_bis and eqref:eq:incoherent_noise, we finally obtain
 #+begin_important
 #+NAME: eq:noise_psd
 \begin{equation}
@@ -4238,6 +4238,254 @@ The generated complementary filters using $\mathcal{H}_\infty$ and the analytica
 - the analytical formula provides a very simple way to generate the complementary filters (and thus the controller), they could even be used to tune the controller online using the parameters $\alpha$ and $\omega_0$. However, these formula have the property that $|H_H|$ and $|H_L|$ are symmetrical with the frequency $\omega_0$ which may not be desirable.
 - while the $\mathcal{H}_\infty$ synthesis of the complementary filters is not as straightforward as using the analytical formula, it provides a more optimized procedure to obtain the complementary filters

-* Bibliography                                                       :ignore:
+* Real World Example of optimal sensor fusion
+** Introduction                                                      :ignore:
+cite:moore19_capac_instr_sensor_fusion_high_bandw_nanop
+
+
+** 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
+
+** Sensor Noise                                                    :noexport:
+#+begin_src matlab
+  A1 = 19.13; % [uV2/Hz]
+  A2 = 0.1632; % [uV2/Hz]
+  A3 = 6.847; % [uV2/Hz]
+  wnc = 3057; % [rad]
+  wx = 7929; % [rad/s]
+ 
+  Fx = 1/(1 - s/wx)/(1 - s/wx);
+  [A B C D] = butter(2, 0.5, 'low');
+  Fx = ss(A, B, C, D);
+
+  Sq = A3*wnc/s + A3;
+  Sx = A1*Fx + A2;
+#+end_src
+
+#+begin_src matlab :exports none
+  freqs = logspace(1, 5, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(Sq, freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(Sx, freqs, 'Hz'))));
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
+#+end_src
+
+** Matlab Code
+Take an Accelerometer and a Geophone both measuring the absolute motion of a structure.
+
+Parameters of the inertial sensors.
+#+begin_src matlab
+  m_acc = 0.01;
+  k_acc = 1e6;
+  c_acc = 20;
+
+  m_geo = 1;
+  k_geo = 1e3;
+  c_geo = 10;
+#+end_src
+
+Transfer function from motion to measurement
+
+For the accelerometer.
+The measurement is the relative motion structure/inertial mass:
+\[ \frac{d}{\ddot{w}} = \frac{-m}{ms^2 + cs + k} \]
+
+For the geophone.
+The measurement is the relative velocity structure/inertial mass:
+\[ \frac{\dot{d}}{\dot{w}} = \frac{-ms^2}{ms^2 + cs + k} \]
+
+#+begin_src matlab
+  G_acc = -m_acc/(m_acc*s^2 + c_acc*s + k_acc); % [m/(m/s^2)]
+  G_geo = -m_geo*s^2/(m_geo*s^2 + c_geo*s + k_geo); % [m/s/m/s]
+#+end_src
+
+Suppose the measure of the relative motion for the accelerometer (capacitive sensor for instance) has a white noise characteristic:
+Suppose the measure of the relative velocity (current flowing through the coil) has a white noise characteristic:
+
+Define the noise characteristics
+#+begin_src matlab
+  n = 1; w0 = 2*pi*5e3; G0 = 5e-12; G1 = 1e-15; Gc = G0/2;
+  L_acc = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
+
+  n = 1; w0 = 2*pi*5e3; G0 = 1e-6; G1 = 1e-8; Gc = G0/2;
+  L_geo = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
+#+end_src
+
+Transfer function of the conversion to obtain the velocity:
+#+begin_src matlab
+  C_acc = (-k_acc/m_acc/(2*pi + s));
+  C_geo = tf(-1);
+#+end_src
+
+Let's plot the noise of both sensors:
+#+begin_src matlab :exports none
+  freqs = logspace(-1, 4, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(L_acc*C_acc, freqs, 'Hz'))), 'DisplayName', 'Acc');
+  plot(freqs, abs(squeeze(freqresp(L_geo*C_geo, freqs, 'Hz'))), 'DisplayName', 'Geo');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Noise ASD [$m/s/\sqrt{Hz}$]');
+  legend('location', 'northeast')
+#+end_src
+
+Dynamics of both sensors
+#+begin_src matlab :exports none
+  freqs = logspace(-1, 4, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(s*G_acc*C_acc, freqs, 'Hz'))), 'DisplayName', 'Acc');
+  plot(freqs, abs(squeeze(freqresp(G_geo*C_geo, freqs, 'Hz'))), 'DisplayName', 'Geo');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude');
+  legend('location', 'northeast')
+#+end_src
+
+** Time domain signals
+#+begin_src matlab
+  Fs = 1e4; % Sampling Frequency [Hz]
+  Ts = 1/Fs; % Sampling Time [s]
+
+  t = 0:Ts:10; % Time Vector [s]
+#+end_src
+
+#+begin_src matlab
+  n_acc = lsim(L_acc*C_acc, sqrt(Fs/2)*randn(length(t), 1), t); % [m/s]
+  n_geo = lsim(L_geo*C_geo, sqrt(Fs/2)*randn(length(t), 1), t); % [m/s]
+#+end_src
+
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(t, n_geo)
+  plot(t, n_acc)
+  hold off;
+#+end_src
+
+** H2 Synthesis
+#+begin_src matlab
+  N1 = L_acc*C_acc;
+  N2 = L_geo*C_geo;
+#+end_src
+
+#+begin_src matlab
+  bodeFig({N1, N2}, logspace(-1, 5, 1000))
+#+end_src
+
+#+begin_src matlab
+  P = [0   N2  1;
+       N1 -N2  0];
+#+end_src
+
+And we do the $\mathcal{H}_2$ synthesis using the =h2syn= command.
+#+begin_src matlab
+  [H1, ~, gamma] = h2syn(P, 1, 1);
+#+end_src
+
+Finally, we define $H_2(s) = 1 - H_1(s)$.
+#+begin_src matlab
+  H2 = 1 - H1;
+#+end_src
+
+#+begin_src matlab
+  bodeFig({H1, H2}, struct('phase', true))
+#+end_src
+
+#+begin_src matlab
+  n_acc_filt = lsim(H1, n_acc, t);
+  n_geo_filt = lsim(H2, n_geo, t);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable([rms(n_acc), rms(n_geo), rms(n_acc_filt + n_geo_filt)]', {'Accelerometer', 'Geophone', 'Super Sensor'}, {'RMS'}, ' %.1e ');
+#+end_src
+
+#+RESULTS:
+|               |     RMS |
+|---------------+---------|
+| Accelerometer | 9.7e-05 |
+| Geophone      | 5.9e-05 |
+| Super Sensor  | 1.5e-05 |
+
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(t, n_geo)
+  plot(t, n_acc)
+  plot(t, n_acc_filt + n_geo_filt)
+  hold off;
+#+end_src
+
+** Signal and Noise
+Velocity Signal:
+#+begin_src matlab
+  v = lsim(1/(1 + s/2/pi/2), 1e-4*sqrt(Fs/2)*randn(length(t), 1), t);
+  v = 1e-4 * sin(2*pi*100*t);
+#+end_src
+
+#+begin_src matlab
+  v_acc = lsim(s*G_acc*C_acc, v, t) + n_acc;
+  v_geo = lsim(G_geo*C_geo,   v, t) + n_geo;
+#+end_src
+
+#+begin_src matlab
+  v_ss = lsim(H1, v_acc, t) + lsim(H2, v_geo, t);
+#+end_src
+
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(t, v_geo)
+  plot(t, v_acc)
+  plot(t, v_ss)
+  plot(t, v, 'k--')
+  hold off;
+  xlim([1, 1+0.1])
+#+end_src
+
+** PSD and CPS
+#+begin_src matlab
+  nx = length(n_acc);
+  na = 16;
+  win = hanning(floor(nx/na));
+
+  [p_acc, f] = pwelch(n_acc, win, 0, [], Fs);
+  [p_geo, ~] = pwelch(n_geo, win, 0, [], Fs);
+  [p_ss,  ~] = pwelch(n_acc_filt + n_geo_filt, win, 0, [], Fs);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(f, p_acc, 'DisplayName', 'Accelerometer');
+  plot(f, p_geo, 'DisplayName', 'Geophone');
+  plot(f, p_ss,  'DisplayName', 'Super Sensor');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]');
+  ylabel('Power Spectral Density $\left[\frac{(m/s)^2}{Hz}\right]$');
+  legend('location', 'southwest');
+#+end_src
+
+** Transfer function of the super sensor
+#+begin_src matlab
+  bodeFig({s*C_acc*G_acc, C_geo*G_geo, s*C_acc*G_acc*H1+C_geo*G_geo*H2}, struct('phase', true))
+#+end_src
+
+* Bibliography                                                        :ignore:
 bibliographystyle:unsrt
 bibliography:ref.bib
--- a/matlab/ref.bib
+++ b/matlab/ref.bib
@@ -181,4 +181,19 @@
  year =         1998,
  doi =          {10.1063/1.1149013},
  url =          {https://doi.org/10.1063/1.1149013},
-}
+}
+
+
+@article{moore19_capac_instr_sensor_fusion_high_bandw_nanop,
+  author          = {Steven Ian Moore and Andrew J. Fleming and Yuen Kuan Yong},
+  title           = {Capacitive Instrumentation and Sensor Fusion for
+                  High-Bandwidth Nanopositioning},
+  journal         = {IEEE Sensors Letters},
+  volume          = 3,
+  number          = 8,
+  pages           = {1-3},
+  year            = 2019,
+  doi             = {10.1109/lsens.2019.2933065},
+  url             = {https://doi.org/10.1109/lsens.2019.2933065},
+  tags            = {complementary filters},
+}