diff --git a/.gitignore b/.gitignore index 8df43c8..f09ce98 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# Emacs +auto/ + # Simulink Real Time *bio.m *pt.m diff --git a/huddle-test-geophones/figs/comp_coherence_formula.png b/huddle-test-geophones/figs/comp_coherence_formula.png new file mode 100644 index 0000000..6e96b60 Binary files /dev/null and b/huddle-test-geophones/figs/comp_coherence_formula.png differ diff --git a/huddle-test-geophones/figs/compare_axis_psd.png b/huddle-test-geophones/figs/compare_axis_psd.png new file mode 100644 index 0000000..94c00a5 Binary files /dev/null and b/huddle-test-geophones/figs/compare_axis_psd.png differ diff --git a/huddle-test-geophones/figs/compare_tf_axis.png b/huddle-test-geophones/figs/compare_tf_axis.png new file mode 100644 index 0000000..921874c Binary files /dev/null and b/huddle-test-geophones/figs/compare_tf_axis.png differ diff --git a/huddle-test-geophones/figs/huddle-test.png b/huddle-test-geophones/figs/huddle-test.png index df2a0bd..bf72425 100644 Binary files a/huddle-test-geophones/figs/huddle-test.png 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 index 62d1797..7c6e819 100644 Binary files a/huddle-test-geophones/figs/huddle_test_results.png and b/huddle-test-geophones/figs/huddle_test_results.png differ diff --git a/huddle-test-geophones/figs/intrumental_noise_V.png b/huddle-test-geophones/figs/intrumental_noise_V.png new file mode 100644 index 0000000..14ba1a8 Binary files /dev/null and b/huddle-test-geophones/figs/intrumental_noise_V.png differ diff --git a/huddle-test-geophones/figs/intrumental_noise_velocity.png b/huddle-test-geophones/figs/intrumental_noise_velocity.png new file mode 100644 index 0000000..c8f2afc Binary files /dev/null and b/huddle-test-geophones/figs/intrumental_noise_velocity.png differ diff --git a/huddle-test-geophones/figs/psd_velocity.png b/huddle-test-geophones/figs/psd_velocity.png index 40716d8..2451e28 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.html b/huddle-test-geophones/index.html index 7d6f41b..311c4a6 100644 Binary files a/huddle-test-geophones/index.html and b/huddle-test-geophones/index.html differ diff --git a/huddle-test-geophones/index.org b/huddle-test-geophones/index.org index 41e3e1f..887b654 100644 --- a/huddle-test-geophones/index.org +++ b/huddle-test-geophones/index.org @@ -43,6 +43,7 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH #+end_src ** Load data +We load the data of the z axis of two geophones. #+begin_src matlab :results none load('mat/data_001.mat', 't', 'x1', 'x2'); dt = t(2) - t(1); @@ -94,7 +95,7 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH #+RESULTS: fig:data_time_domain_zoom [[file:figs/data_time_domain_zoom.png]] -** Compute PSD +** Computation of the ASD of the measured voltage We first define the parameters for the frequency domain analysis. #+begin_src matlab :results none win = hanning(ceil(length(x1)/100)); @@ -106,15 +107,17 @@ We first define the parameters for the frequency domain analysis. [pxx2, ~] = pwelch(x2, win, [], [], Fs); #+end_src -** Take into account sensibility of Geophone +** Scaling to take into account the sensibility of the geophone and the voltage amplifier The Geophone used are L22. +Their sensibility are shown on figure [[fig:geophone_sensibility]]. + #+begin_src matlab :results none S0 = 88; % Sensitivity [V/(m/s)] f0 = 2; % Cut-off frequnecy [Hz] S = (s/2/pi/f0)/(1+s/2/pi/f0); #+end_src -#+begin_src matlab :results none +#+begin_src matlab :results none :exports none figure; bodeFig({S}); ylabel('Amplitude [V/(m/s)]') @@ -132,11 +135,8 @@ The Geophone used are L22. [[file:figs/geophone_sensibility.png]] -We take into account the gain of the electronics. -The cut-off frequency is set at 1kHz. - -- [ ] Check what is the order of the filter -- [ ] Maybe I should not use this filter as there is no high frequencies anyway? +We also take into account the gain of the electronics which is here set to be $60dB$. +The amplifiers also include a low pass filter with a cut-off frequency set at 1kHz. #+begin_src matlab :results none G0 = 60; % [dB] @@ -144,11 +144,21 @@ The cut-off frequency is set at 1kHz. G = G0/(1+s/2/pi/1000); #+end_src +We divide the ASD measured (in $\text{V}/\sqrt{\text{Hz}}$) by the transfer function of the voltage amplifier to obtain the ASD of the voltage across the geophone. +We further divide the result by the sensibility of the Geophone to obtain the ASD of the velocity in $m/s/\sqrt{Hz}$. + +#+begin_src matlab :results none + scaling = 1./squeeze(abs(freqresp(G, f, 'Hz')))./squeeze(abs(freqresp(S, f, 'Hz'))); +#+end_src + +** Computation of the ASD of the velocity +The ASD of the measured velocity is shown on figure [[fig:psd_velocity]]. + #+begin_src matlab :results none figure; hold on; - 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')))); + plot(f, sqrt(pxx1)./scaling); + plot(f, sqrt(pxx2)./scaling); hold off; set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); @@ -168,11 +178,16 @@ The cut-off frequency is set at 1kHz. [[file:figs/psd_velocity.png]] ** Transfer function between the two geophones +We here compute the transfer function from one geophone to the other. +The result is shown on figure [[fig:tf_geophones]]. + +We also compute the coherence between the two signals (figure [[fig:coh_geophones]]). + #+begin_src matlab :results none [T12, ~] = tfestimate(x1, x2, win, [], [], Fs); #+end_src -#+begin_src matlab :results none +#+begin_src matlab :results none :exports none figure; ax1 = subplot(2, 1, 1); plot(f, abs(T12)); @@ -206,7 +221,7 @@ The cut-off frequency is set at 1kHz. [coh12, ~] = mscohere(x1, x2, win, [], [], Fs); #+end_src -#+begin_src matlab :results none +#+begin_src matlab :results none :exports none figure; plot(f, coh12); set(gca, 'xscale', 'log'); @@ -225,63 +240,62 @@ The cut-off frequency is set at 1kHz. #+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) {}; +** Estimation of the sensor noise +The technique to estimate the sensor noise is taken from cite:barzilai98_techn_measur_noise_sensor_presen. - \draw[] (U) node[above right]{$U$} -- ++(1, 0) node[]{$\bullet$}; - \draw[->] ($(U)+(1, 0)$) |- (S1.west); - \draw[->] ($(U)+(1, 0)$) |- (S2.west); +The coherence between signals $X$ and $Y$ is defined as follow +\[ \gamma^2_{XY}(\omega) = \frac{|G_{XY}(\omega)|^2}{|G_{X}(\omega)| |G_{Y}(\omega)|} \] +where $|G_X(\omega)|$ is the output Power Spectral Density (PSD) of signal $X$ and $|G_{XY}(\omega)|$ is the Cross Spectral Density (CSD) of signal $X$ and $Y$. - \draw[->] (S1.east) -- (add1.west); - \draw[->] (S2.east) -- (add2.west); +The PSD and CSD are defined as follow: +\begin{align} + |G_X(\omega)| &= \frac{2}{n_d T} \sum^{n_d}_{n=1} \left| X_k(\omega, T) \right|^2 \\ + |G_{XY}(\omega)| &= \frac{2}{n_d T} \sum^{n_d}_{n=1} [ X_k^*(\omega, T) ] [ Y_k(\omega, T) ] +\end{align} +where: +- $n_d$ is the number for records averaged +- $T$ is the length of each record +- $X_k(\omega, T)$ is the finite Fourier transform of the kth record +- $X_k^*(\omega, T)$ is its complex conjugate - \draw[->] (add1.east) -- ++(1, 0) node[above]{$X_1$}; - \draw[->] (add2.east) -- ++(1, 0) node[above]{$X_2$}; +The =mscohere= function is compared with this formula on Appendix (section [[sec:coherence]]), it is shown that it is identical. - \draw[<-] (add1.north) -- ++(0, 0.8)node[right]{$N_1$}; - \draw[<-] (add2.north) -- ++(0, 0.8)node[right]{$N_2$}; - \end{tikzpicture} -#+end_src +Figure [[fig:huddle_test]] illustrate a block diagram model of the system used to determine the sensor noise of the geophone. + +Two geophones are mounted side by side to ensure that they are exposed by the same motion input $U$. + +Each sensor has noise $N$ and $M$. #+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$. +We here assume that each sensor has the same magnitude of instrumental noise ($N = M$). +We also assume that $H_1 = H_2 = 1$. -\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) \] +We then obtain: +#+NAME: eq:coh_bis +\begin{equation} + \gamma_{XY}^2(\omega) = \frac{1}{1 + 2 \left( \frac{|G_N(\omega)|}{|G_U(\omega)|} \right) + \left( \frac{|G_N(\omega)|}{|G_U(\omega)|} \right)^2} +\end{equation} +Since the input signal $U$ and the instrumental noise $N$ are incoherent: +#+NAME: eq:incoherent_noise +\begin{equation} + |G_X(\omega)| = |G_N(\omega)| + |G_U(\omega)| +\end{equation} +From equations [[eq:coh_bis]] and [[eq:incoherent_noise]], we finally obtain +#+begin_important +#+NAME: eq:noise_psd +\begin{equation} + |G_N(\omega)| = |G_X(\omega)| \left( 1 - \sqrt{\gamma_{XY}^2(\omega)} \right) +\end{equation} +#+end_important +The instrumental noise is computed below. The result in V^2/Hz is shown on figure [[fig:intrumental_noise_V]]. #+begin_src matlab :results none - S = abs(T12.*pxx1); - - N = pxx2 - (T12.^2).*pxx1; - N = abs(N)/2; + pxxN = pxx1.*(1 - coh12); #+end_src #+begin_src matlab :results none @@ -289,20 +303,184 @@ and hold on; plot(f, pxx1, '-'); plot(f, pxx2, '-'); - plot(f, N, 'k:', 'linewidth', 1); + plot(f, pxxN, 'k--'); hold off; set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('PSD [$V^2/Hz$]'); xlim([1, 500]); - legend('$\Phi_{ss} (f)$','$\Phi_{nn} (f)$') #+end_src -#+NAME: fig:huddle_test_results +#+NAME: fig:intrumental_noise_V #+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") +#+begin_src matlab :var filepath="figs/intrumental_noise_V.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]] +#+NAME: fig:intrumental_noise_V +#+CAPTION: Instrumental Noise and Measurement in $V^2/Hz$ +#+RESULTS: fig:intrumental_noise_V +[[file:figs/intrumental_noise_V.png]] + +This is then further converted into velocity and compared with the ground velocity measurement. (figure [[fig:intrumental_noise_velocity]]) +#+begin_src matlab :results none + figure; + hold on; + plot(f, sqrt(pxx1).*scaling, '-'); + plot(f, sqrt(pxx2).*scaling, '-'); + plot(f, sqrt(pxxN).*scaling, 'k--'); + hold off; + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('PSD [$m/s/\sqrt{Hz}$]'); + xlim([1, 500]); +#+end_src + +#+NAME: fig:intrumental_noise_velocity +#+HEADER: :tangle no :exports results :results value raw replace :noweb yes +#+begin_src matlab :var filepath="figs/intrumental_noise_velocity.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:intrumental_noise_velocity +#+CAPTION: Instrumental Noise and Measurement in $m/s/\sqrt{Hz}$ +#+RESULTS: fig:intrumental_noise_velocity +[[file:figs/intrumental_noise_velocity.png]] + +* Compare axis +** Matlab Init :noexport:ignore: +#+begin_src matlab :exports none :results silent :noweb yes + <> +#+end_src + +** Load data +#+begin_src matlab :results none + z = load('mat/data_001.mat', 't', 'x1', 'x2'); + east = load('mat/data_002.mat', 't', 'x1', 'x2'); + north = load('mat/data_003.mat', 't', 'x1', 'x2'); +#+end_src + +** Compare PSD +#+begin_src matlab :results none + [pz1, fz] = pwelch(z.x1, hanning(ceil(length(z.x1)/100)), [], [], 1/(z.t(2)-z.t(1))); + [pz2, ~] = pwelch(z.x2, hanning(ceil(length(z.x2)/100)), [], [], 1/(z.t(2)-z.t(1))); + + [pe1, fe] = pwelch(east.x1, hanning(ceil(length(east.x1)/100)), [], [], 1/(east.t(2)-east.t(1))); + [pe2, ~] = pwelch(east.x2, hanning(ceil(length(east.x2)/100)), [], [], 1/(east.t(2)-east.t(1))); + + [pn1, fn] = pwelch(north.x1, hanning(ceil(length(north.x1)/100)), [], [], 1/(north.t(2)-north.t(1))); + [pn2, ~] = pwelch(north.x2, hanning(ceil(length(north.x2)/100)), [], [], 1/(north.t(2)-north.t(1))); +#+end_src + +#+begin_src matlab :results none :exports none + figure; + hold on; + plot(fz, sqrt(pz1), '-', 'Color', [0 0.4470 0.7410], 'DisplayName', 'z'); + plot(fz, sqrt(pz2), '--', 'Color', [0 0.4470 0.7410], 'HandleVisibility', 'off'); + plot(fe, sqrt(pe1), '-', 'Color', [0.8500 0.3250 0.0980], 'DisplayName', 'east'); + plot(fe, sqrt(pe2), '--', 'Color', [0.8500 0.3250 0.0980], 'HandleVisibility', 'off'); + plot(fn, sqrt(pn1), '-', 'Color', [0.9290 0.6940 0.1250], 'DisplayName', 'north'); + plot(fn, sqrt(pn2), '--', 'Color', [0.9290 0.6940 0.1250], 'HandleVisibility', 'off'); + hold off; + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]'); + legend('Location', 'northeast'); + xlim([0, 500]); +#+end_src + +#+NAME: fig:compare_axis_psd +#+HEADER: :tangle no :exports results :results value raw replace :noweb yes +#+begin_src matlab :var filepath="figs/compare_axis_psd.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:compare_axis_psd +#+CAPTION: Compare the measure PSD of the two geophones for the three axis +#+RESULTS: fig:compare_axis_psd +[[file:figs/compare_axis_psd.png]] + +** Compare TF +#+begin_src matlab :results none + [Tz, fz] = tfestimate(z.x1, z.x2, hanning(ceil(length(z.x1)/100)), [], [], 1/(z.t(2)-z.t(1))); + [Te, fe] = tfestimate(east.x1, east.x2, hanning(ceil(length(east.x1)/100)), [], [], 1/(east.t(2)-east.t(1))); + [Tn, fn] = tfestimate(north.x1, north.x2, hanning(ceil(length(north.x1)/100)), [], [], 1/(north.t(2)-north.t(1))); +#+end_src + +#+begin_src matlab :results none :exports none + figure; + ax1 = subplot(2, 1, 1); + hold on; + plot(fz, abs(Tz), 'DisplayName', 'z'); + plot(fe, abs(Te), 'DisplayName', 'east'); + plot(fn, abs(Tn), 'DisplayName', 'north'); + hold off; + set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log'); + set(gca, 'XTickLabel',[]); + ylabel('Magnitude'); + legend('Location', 'southwest'); + + ax2 = subplot(2, 1, 2); + hold on; + plot(fz, mod(180+180/pi*phase(Tz), 360)-180); + plot(fe, mod(180+180/pi*phase(Te), 360)-180); + plot(fn, mod(180+180/pi*phase(Tn), 360)-180); + hold off; + set(gca, 'xscale', 'log'); + ylim([-180, 180]); + yticks([-180, -90, 0, 90, 180]); + xlabel('Frequency [Hz]'); ylabel('Phase'); + + linkaxes([ax1,ax2],'x'); + xlim([1, 500]); +#+end_src + +#+NAME: fig:compare_tf_axis +#+HEADER: :tangle no :exports results :results value raw replace :noweb yes +#+begin_src matlab :var filepath="figs/compare_tf_axis.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:compare_tf_axis +#+CAPTION: Compare the transfer function from one geophone to the other for the 3 axis +#+RESULTS: fig:compare_tf_axis +[[file:figs/compare_tf_axis.png]] +* Appendix +** Computation of coherence from PSD and CSD + <> +#+begin_src matlab :results none + load('mat/data_001.mat', 't', 'x1', 'x2'); + dt = t(2) - t(1); + Fs = 1/dt; + win = hanning(ceil(length(x1)/100)); +#+end_src + +#+begin_src matlab :results none + pxy = cpsd(x1, x2, win, [], [], Fs); + pxx = pwelch(x1, win, [], [], Fs); + pyy = pwelch(x2, win, [], [], Fs); + coh = mscohere(x1, x2, win, [], [], Fs); +#+end_src + +#+begin_src matlab :results none + figure; + hold on; + plot(f, abs(pxy).^2./abs(pxx)./abs(pyy), '-'); + plot(f, coh, '--'); + hold off; + set(gca, 'xscale', 'log'); + xlabel('Frequency'); ylabel('Coherence'); + xlim([1, 500]); +#+end_src + +#+NAME: fig:comp_coherence_formula +#+HEADER: :tangle no :exports results :results value raw replace :noweb yes +#+begin_src matlab :var filepath="figs/comp_coherence_formula.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png") + <> +#+end_src + +#+NAME: fig:comp_coherence_formula +#+CAPTION: Comparison of =mscohere= and manual computation +#+RESULTS: fig:comp_coherence_formula +[[file:figs/comp_coherence_formula.png]] + +* Bibliography :ignore: +bibliographystyle:unsrt +bibliography:ref.bib diff --git a/huddle-test-geophones/readme.org b/huddle-test-geophones/readme.org index e478308..aab2d19 100644 --- a/huddle-test-geophones/readme.org +++ b/huddle-test-geophones/readme.org @@ -1,8 +1,9 @@ -* DONE Register data on the computer - CLOSED: [2019-04-17 mer. 17:26] +* TODO [#B] Find the documentation of the amplifier to know the order of the filters +* TODO [#A] Shake a little bit the geophones to see if we have better measurements on X and Y axis +* Measurements -* Measurements: - -| data_001.mat | Z axis | -| data_002.mat | East | -| data_003.mat | North | \ No newline at end of file +| Filename | Description | +|--------------+-------------| +| data_001.mat | Z axis | +| data_002.mat | East | +| data_003.mat | North | diff --git a/huddle-test-geophones/ref.bib b/huddle-test-geophones/ref.bib new file mode 100644 index 0000000..a602ec3 --- /dev/null +++ b/huddle-test-geophones/ref.bib @@ -0,0 +1,29 @@ +@article{barzilai98_techn_measur_noise_sensor_presen, + author = {Aaron Barzilai and Tom VanZandt and Tom Kenny}, + title = {Technique for Measurement of the Noise of a Sensor in the + Presence of Large Background Signals}, + journal = {Review of Scientific Instruments}, + volume = 69, + number = 7, + pages = {2767-2772}, + year = 1998, + doi = {10.1063/1.1149013}, + url = {https://doi.org/10.1063/1.1149013}, +} + + +@article{kirchhoff17_huddl_test_measur_near_johns, + author = {R. Kirchhoff and C. M. Mow-Lowry and V. B. Adya and G. + Bergmann and S. Cooper and M. M. Hanke and P. Koch and S. M. + K{\"o}hlenbeck and J. Lehmann and P. Oppermann and J. + W{\"o}hler and D. S. Wu and H. L{\"u}ck and K. A. Strain}, + title = {Huddle Test Measurement of a Near Johnson Noise Limited + Geophone}, + journal = {Review of Scientific Instruments}, + volume = 88, + number = 11, + pages = 115008, + year = 2017, + doi = {10.1063/1.5000592}, + url = {https://doi.org/10.1063/1.5000592}, +} \ No newline at end of file