nass-micro-station-measurem.../huddle-test-geophones/index.org

487 lines
16 KiB
Org Mode

#+TITLE:SpeedGoat
:DRAWER:
#+STARTUP: overview
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="../js/readtheorg.js"></script>
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :results output
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:
* Setup
Two L22 geophones are used.
They are placed on the ID31 granite.
They are leveled.
The signals are amplified using voltage amplifier with a gain of 60dB.
The voltage amplifiers include a low pass filter with a cut-off frequency at 1kHz.
#+name: fig:figure_name
#+caption: Setup
#+attr_html: :width 500px
[[file:./figs/setup.jpg]]
#+name: fig:figure_name
#+caption: Geophones
#+attr_html: :width 500px
[[file:./figs/geophones.jpg]]
* Signal Processing
** Matlab Init :noexport:ignore:
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+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);
#+end_src
** Time Domain Data
#+begin_src matlab :results none
figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([t(1), t(end)]);
#+end_src
#+NAME: fig:data_time_domain
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/data_time_domain.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:data_time_domain
#+CAPTION: Time domain Data
#+RESULTS: fig:data_time_domain
[[file:figs/data_time_domain.png]]
#+begin_src matlab :results none
figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([0 1]);
#+end_src
#+NAME: fig:data_time_domain_zoom
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/data_time_domain_zoom.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:data_time_domain_zoom
#+CAPTION: Time domain Data - Zoom
#+RESULTS: fig:data_time_domain_zoom
[[file:figs/data_time_domain_zoom.png]]
** 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));
Fs = 1/dt;
#+end_src
#+begin_src matlab :results none
[pxx1, f] = pwelch(x1, win, [], [], Fs);
[pxx2, ~] = pwelch(x2, win, [], [], Fs);
#+end_src
** 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 :exports none
figure;
bodeFig({S});
ylabel('Amplitude [V/(m/s)]')
#+end_src
#+NAME: fig:geophone_sensibility
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/geophone_sensibility.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:geophone_sensibility
#+CAPTION: Sensibility of the Geophone
#+RESULTS: fig:geophone_sensibility
[[file:figs/geophone_sensibility.png]]
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]
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)./scaling);
plot(f, sqrt(pxx2)./scaling);
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
xlim([2, 500]);
#+end_src
#+NAME: fig:psd_velocity
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/psd_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:psd_velocity
#+CAPTION: Spectral density of the velocity
#+RESULTS: fig:psd_velocity
[[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 :exports none
figure;
ax1 = subplot(2, 1, 1);
plot(f, abs(T12));
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude');
ax2 = subplot(2, 1, 2);
plot(f, mod(180+180/pi*phase(T12), 360)-180);
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:tf_geophones
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/tf_geophones.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+NAME: fig:tf_geophones
#+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 :exports 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]]
** Estimation of the sensor noise
The technique to estimate the sensor noise is taken from cite:barzilai98_techn_measur_noise_sensor_presen.
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$.
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
The =mscohere= function is compared with this formula on Appendix (section [[sec:coherence]]), it is shown that it is identical.
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
[[file:figs/huddle-test.png]]
We here assume that each sensor has the same magnitude of instrumental noise ($N = M$).
We also assume that $H_1 = H_2 = 1$.
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
pxxN = pxx1.*(1 - coh12);
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, pxx1, '-');
plot(f, pxx2, '-');
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]);
#+end_src
#+NAME: fig:intrumental_noise_V
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/intrumental_noise_V.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
<<plt-matlab>>
#+end_src
#+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")
<<plt-matlab>>
#+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
<<matlab-init>>
#+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")
<<plt-matlab>>
#+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")
<<plt-matlab>>
#+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
<<sec:coherence>>
#+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")
<<plt-matlab>>
#+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