Add huddle test

huddle-test-geophones/index.org

@@ -95,9 +95,15 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH
 [[file:figs/data_time_domain_zoom.png]]
 
 ** Compute PSD
+We first define the parameters for the frequency domain analysis.
 #+begin_src matlab :results none
-  [pxx1, f1] = pwelch(x1, hanning(ceil(1/dt)), 0, [], 1/dt);
-  [pxx2, f2] = pwelch(x2, hanning(ceil(1/dt)), 0, [], 1/dt);
+  win = hanning(ceil(length(x1)/100));
+  Fs = 1/dt;
+#+end_src
+
+#+begin_src matlab :results none
+  [pxx1, f] = pwelch(x1, win, [], [], Fs);
+  [pxx2, ~] = pwelch(x2, win, [], [], Fs);
 #+end_src
 
 ** Take into account sensibility of Geophone
@@ -141,8 +147,8 @@ The cut-off frequency is set at 1kHz.
 #+begin_src matlab :results none
   figure;
   hold on;
-  plot(f1, sqrt(pxx1)./squeeze(abs(freqresp(G, f1, 'Hz')))./squeeze(abs(freqresp(S, f1, 'Hz'))));
-  plot(f2, sqrt(pxx2)./squeeze(abs(freqresp(G, f2, 'Hz')))./squeeze(abs(freqresp(S, f2, 'Hz'))));
+  plot(f, sqrt(pxx1)./squeeze(abs(freqresp(G, f, 'Hz')))./squeeze(abs(freqresp(S, f1, 'Hz'))));
+  plot(f, sqrt(pxx2)./squeeze(abs(freqresp(G, f, 'Hz')))./squeeze(abs(freqresp(S, f2, 'Hz'))));
   hold off;
   set(gca, 'xscale', 'log');
   set(gca, 'yscale', 'log');
@@ -160,21 +166,22 @@ The cut-off frequency is set at 1kHz.
 #+CAPTION: Spectral density of the velocity
 #+RESULTS: fig:psd_velocity
 [[file:figs/psd_velocity.png]]
+
 ** Transfer function between the two geophones
 #+begin_src matlab :results none
-  [T12, f12] = tfestimate(x1, x2, hanning(ceil(length(x1)/100)), [], [], 1/dt);
+  [T12, ~] = tfestimate(x1, x2, win, [], [], Fs);
 #+end_src
 
 #+begin_src matlab :results none
   figure;
   ax1 = subplot(2, 1, 1);
-  plot(f12, abs(T12));
+  plot(f, abs(T12));
   set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
   set(gca, 'XTickLabel',[]);
   ylabel('Magnitude');
 
   ax2 = subplot(2, 1, 2);
-  plot(f12, mod(180+180/pi*phase(T12), 360)-180);
+  plot(f, mod(180+180/pi*phase(T12), 360)-180);
   set(gca, 'xscale', 'log');
   ylim([-180, 180]);
   yticks([-180, -90, 0, 90, 180]);
@@ -194,3 +201,108 @@ The cut-off frequency is set at 1kHz.
 #+CAPTION: Estimated transfer function between the two geophones
 #+RESULTS: fig:tf_geophones
 [[file:figs/tf_geophones.png]]
+
+#+begin_src matlab :results none
+  [coh12, ~] = mscohere(x1, x2, win, [], [], Fs);
+#+end_src
+
+#+begin_src matlab :results none
+  figure;
+  plot(f, coh12);
+  set(gca, 'xscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Coherence');
+  ylim([0,1]); xlim([1, 500]);
+#+end_src
+
+#+NAME: fig:coh_geophones
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/coh_geophones.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:coh_geophones
+#+CAPTION: Cohererence between the signals of the two geophones
+#+RESULTS: fig:coh_geophones
+[[file:figs/coh_geophones.png]]
+
+** Huddle Test
+#+NAME: fig:huddle_test
+#+HEADER: :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
+#+HEADER: :imagemagick t :fit yes :iminoptions -scale 100% -density 150 :imoutoptions -quality 100
+#+HEADER: :results raw replace :buffer no :eval no-export :exports both :mkdirp yes
+#+HEADER: :output-dir figs
+#+begin_src latex :file huddle-test.pdf :post pdf2svg(file=*this*, ext="png") :exports results
+  \begin{tikzpicture}
+    \coordinate[] (U) at (0, 0) {};
+    \node[block, above right=0.5 and 2 of U] (S1) {$S_1$};
+    \node[block, below right=0.5 and 2 of U] (S2) {$S_2$};
+    \node[addb={+}{}{}{}{}, right=0.5 of S1] (add1) {};
+    \node[addb={+}{}{}{}{}, right=0.5 of S2] (add2) {};
+
+    \draw[] (U) node[above right]{$U$} -- ++(1, 0) node[]{$\bullet$};
+    \draw[->] ($(U)+(1, 0)$)  |- (S1.west);
+    \draw[->] ($(U)+(1, 0)$)  |- (S2.west);
+
+    \draw[->] (S1.east)  -- (add1.west);
+    \draw[->] (S2.east)  -- (add2.west);
+
+    \draw[->] (add1.east)  -- ++(1, 0) node[above]{$X_1$};
+    \draw[->] (add2.east)  -- ++(1, 0) node[above]{$X_2$};
+
+    \draw[<-] (add1.north)  -- ++(0, 0.8)node[right]{$N_1$};
+    \draw[<-] (add2.north)  -- ++(0, 0.8)node[right]{$N_2$};
+  \end{tikzpicture}
+#+end_src
+
+#+NAME: fig:huddle_test
+#+CAPTION: Huddle test block diagram
+#+RESULTS: fig:huddle_test
+[[file:figs/huddle-test.png]]
+
+We are measuring $X_1$ and $X_2$.
+The goal is to determine $N$.
+
+\begin{align*}
+  X_1(\omega) &= S_1(\omega) U(\omega) + N_1(\omega)\\
+  X_2(\omega) &= S_2(\omega) U(\omega) + N_2(\omega)
+\end{align*}
+
+Then
+\[ X_2(\omega) = \frac{S_2(\omega)}{S_1(\omega)} X_1(\omega) + N_2(\omega) - \frac{S_2(\omega)}{S_1(\omega)}N_1(\omega) \]
+
+We suppose $N_1 = N_2 = N$
+\[ N_2(\omega) - \frac{S_2(\omega)}{S_1(\omega)}N_1(\omega) = \left( 1 - \frac{S_2(\omega)}{S_1(\omega)}\right) N(\omega) \]
+and
+\[ N(\omega) = \frac{S_2(\omega)}{S_1(\omega)} X_2(\omega) - N_2(\omega) - \frac{S_2(\omega)}{S_1(\omega)}N_1(\omega) \]
+
+
+
+#+begin_src matlab :results none
+  S = abs(T12.*pxx1);
+
+  N = pxx2 - (T12.^2).*pxx1;
+  N = abs(N)/2;
+#+end_src
+
+#+begin_src matlab :results none
+  figure;
+  hold on;
+  plot(f, pxx1, '-');
+  plot(f, pxx2, '-');
+  plot(f, N, 'k:', 'linewidth', 1);
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlim([1, 500]);
+  legend('$\Phi_{ss} (f)$','$\Phi_{nn} (f)$')
+#+end_src
+
+#+NAME: fig:huddle_test_results
+#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
+#+begin_src matlab :var filepath="figs/huddle_test_results.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
+  <<plt-matlab>>
+#+end_src
+
+#+NAME: fig:huddle_test_results
+#+CAPTION: Results of the Huddle test
+#+RESULTS: fig:huddle_test_results
+[[file:figs/huddle_test_results.png]]