diff --git a/huddle-test-geophones/figs/coh_geophones.png b/huddle-test-geophones/figs/coh_geophones.png new file mode 100644 index 0000000..ce75596 Binary files /dev/null and b/huddle-test-geophones/figs/coh_geophones.png differ diff --git a/huddle-test-geophones/figs/huddle-test.png b/huddle-test-geophones/figs/huddle-test.png new file mode 100644 index 0000000..df2a0bd Binary files /dev/null and b/huddle-test-geophones/figs/huddle-test.png differ diff --git a/huddle-test-geophones/figs/huddle_test_results.png b/huddle-test-geophones/figs/huddle_test_results.png new file mode 100644 index 0000000..62d1797 Binary files /dev/null and b/huddle-test-geophones/figs/huddle_test_results.png differ diff --git a/huddle-test-geophones/figs/psd_velocity.png b/huddle-test-geophones/figs/psd_velocity.png index a3ec32d..40716d8 100644 Binary files a/huddle-test-geophones/figs/psd_velocity.png and b/huddle-test-geophones/figs/psd_velocity.png differ diff --git a/huddle-test-geophones/index.org b/huddle-test-geophones/index.org index b6da652..41e3e1f 100644 --- a/huddle-test-geophones/index.org +++ b/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") + <> +#+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") + <> +#+end_src + +#+NAME: fig:huddle_test_results +#+CAPTION: Results of the Huddle test +#+RESULTS: fig:huddle_test_results +[[file:figs/huddle_test_results.png]]