609 lines
20 KiB
Org Mode
609 lines
20 KiB
Org Mode
#+TITLE:Huddle Test of the L22 Geophones
|
|
:DRAWER:
|
|
#+STARTUP: overview
|
|
|
|
#+LANGUAGE: en
|
|
#+EMAIL: dehaeze.thomas@gmail.com
|
|
#+AUTHOR: Dehaeze Thomas
|
|
|
|
#+HTML_LINK_HOME: ../index.html
|
|
#+HTML_LINK_UP: ../index.html
|
|
|
|
#+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>
|
|
|
|
#+HTML_MATHJAX: align: center tagside: right font: TeX
|
|
|
|
#+PROPERTY: header-args:matlab :session *MATLAB*
|
|
#+PROPERTY: header-args:matlab+ :comments org
|
|
#+PROPERTY: header-args:matlab+ :results none
|
|
#+PROPERTY: header-args:matlab+ :exports both
|
|
#+PROPERTY: header-args:matlab+ :eval no-export
|
|
#+PROPERTY: header-args:matlab+ :output-dir figs
|
|
|
|
#+PROPERTY: header-args:shell :eval no-export
|
|
:END:
|
|
|
|
* Experimental 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 includes:
|
|
- an high pass filter with a cut-off frequency at 1.5Hz (AC option)
|
|
- a low pass filter with a cut-off frequency at 1kHz
|
|
|
|
#+name: fig:figure_name
|
|
#+caption: Setup
|
|
#+attr_html: :width 500px
|
|
[[file:./img/setup.jpg]]
|
|
|
|
#+name: fig:figure_name
|
|
#+caption: Geophones
|
|
#+attr_html: :width 500px
|
|
[[file:./img/geophones.jpg]]
|
|
|
|
* Signal Processing
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/huddle_test_signal_processing.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:huddle_test_signal_processing>>
|
|
|
|
** ZIP file containing the data and matlab files :ignore:
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/huddle_test_signal_processing.m -nt data/huddle_test_signal_processing.zip ]; then
|
|
cp matlab/huddle_test_signal_processing.m huddle_test_signal_processing.m;
|
|
zip data/huddle_test_signal_processing \
|
|
mat/data_001.mat \
|
|
huddle_test_signal_processing.m;
|
|
rm huddle_test_signal_processing.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/huddle_test_signal_processing.zip][here]].
|
|
#+end_note
|
|
|
|
** 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
|
|
|
|
** 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
|
|
Fs = 1/dt; % [Hz]
|
|
|
|
win = hanning(ceil(10*Fs));
|
|
#+end_src
|
|
|
|
Then we compute the Power Spectral Density using =pwelch= function.
|
|
#+begin_src matlab :results none
|
|
[pxx1, f] = pwelch(x1, win, [], [], Fs);
|
|
[pxx2, ~] = pwelch(x2, win, [], [], Fs);
|
|
#+end_src
|
|
|
|
And we plot the result on figure [[fig:asd_voltage]].
|
|
|
|
#+begin_src matlab :results none
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx1));
|
|
plot(f, sqrt(pxx2));
|
|
hold off;
|
|
set(gca, 'xscale', 'log');
|
|
set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD of the measured Voltage $\left[\frac{V}{\sqrt{Hz}}\right]$')
|
|
xlim([0.1, 500]);
|
|
#+end_src
|
|
|
|
#+NAME: fig:asd_voltage
|
|
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
|
|
#+begin_src matlab :var filepath="figs/asd_voltage.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:asd_voltage
|
|
#+CAPTION: Amplitude Spectral Density of the measured voltage
|
|
#+RESULTS: fig:asd_voltage
|
|
[[file:figs/asd_voltage.png]]
|
|
|
|
** Scaling to take into account the sensibility of the geophone and the voltage amplifier
|
|
The Geophone used are L22. Their sensibility is shown on figure [[fig:geophone_sensibility]].
|
|
|
|
#+begin_src matlab :results none
|
|
S0 = 88; % Sensitivity [V/(m/s)]
|
|
f0 = 2; % Cut-off frequnecy [Hz]
|
|
|
|
S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :results none :exports none
|
|
figure;
|
|
bodeFig({S}, logspace(-1, 2, 1000));
|
|
ylabel('Amplitude $\left[\frac{V}{m/s}\right]$')
|
|
#+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$.
|
|
|
|
#+begin_src matlab :results none
|
|
G0_db = 60; % [dB]
|
|
|
|
G0 = 10^(G0_db/20); % [abs]
|
|
#+end_src
|
|
|
|
We divide the ASD measured (in $\text{V}/\sqrt{\text{Hz}}$) by the gain 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(G0*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('ASD of the measured Velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
|
|
xlim([0.1, 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: Amplitude Spectral Density of the Velocity
|
|
#+RESULTS: fig:psd_velocity
|
|
[[file:figs/psd_velocity.png]]
|
|
|
|
We also plot the ASD in displacement (figure [[fig:asd_displacement]]);
|
|
|
|
#+begin_src matlab :results none
|
|
figure;
|
|
hold on;
|
|
plot(f, (sqrt(pxx1).*scaling)./(2*pi*f));
|
|
plot(f, (sqrt(pxx2).*scaling)./(2*pi*f));
|
|
hold off;
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
|
|
xlim([0.1, 500]);
|
|
#+end_src
|
|
|
|
#+NAME: fig:asd_displacement
|
|
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
|
|
#+begin_src matlab :var filepath="figs/asd_displacement.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
|
|
<<plt-matlab>>
|
|
#+end_src
|
|
|
|
#+NAME: fig:asd_displacement
|
|
#+CAPTION: Amplitude Spectral Density of the Displacement
|
|
#+RESULTS: fig:asd_displacement
|
|
[[file:figs/asd_displacement.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 [deg]');
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([0.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([0.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 $S_1 = S_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 of the measured Voltage $\left[\frac{V^2}{Hz}\right]$');
|
|
xlim([0.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="full-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('ASD of the Velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$');
|
|
xlim([0.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="full-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
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/huddle_test_compare_axis.m
|
|
:header-args:matlab+: :comments org :mkdirp yes
|
|
:END:
|
|
<<sec:huddle_test_compare_axis>>
|
|
|
|
#+begin_src bash :exports none :results none
|
|
if [ matlab/huddle_test_compare_axis.m -nt data/huddle_test_compare_axis.zip ]; then
|
|
cp matlab/huddle_test_compare_axis.m huddle_test_compare_axis.m;
|
|
zip data/huddle_test_compare_axis \
|
|
mat/data_001.mat \
|
|
mat/data_002.mat \
|
|
mat/data_003.mat \
|
|
huddle_test_compare_axis.m;
|
|
rm huddle_test_compare_axis.m;
|
|
fi
|
|
#+end_src
|
|
|
|
#+begin_note
|
|
All the files (data and Matlab scripts) are accessible [[file:data/huddle_test_compare_axis.zip][here]].
|
|
#+end_note
|
|
|
|
** 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
|
|
|
|
** Load data
|
|
We first load the data for the three axis.
|
|
#+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
|
|
The PSD for each axis of the two geophones are computed.
|
|
#+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
|
|
|
|
We compare them. The result is shown on figure [[fig:compare_axis_psd]].
|
|
#+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
|
|
The transfer functions from one geophone to the other are also computed for each axis.
|
|
The result is shown on figure [[fig:compare_tf_axis]].
|
|
|
|
#+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
|