tdehaeze/nass-micro-station-measurements
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Add links to matlab and mat files

Browse Source
This commit is contained in:
Thomas Dehaeze 2019-04-18 17:11:25 +02:00
parent 97755ac9f0
commit 36f1e7d875
4 changed files with 446 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

81
huddle-test-geophones/compare_axis.m Normal file
View File

@ -0,0 +1,81 @@
% Matlab Init :noexport:ignore:
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Initialize ans with org-babel
ans = 0;
% Load data
% We first load the data for the three axis.
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');
% Compare PSD
% The PSD for each axis of the two geophones are computed.
[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)));
% We compare them. The result is shown on figure [[fig:compare_axis_psd]].
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]);
% 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]].
[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)));
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]);

179
huddle-test-geophones/index.html
View File

@ -3,7 +3,7 @@
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2019-04-18 jeu. 17:02 -->
<!-- 2019-04-18 jeu. 17:11 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>SpeedGoat</title>
@ -276,36 +276,36 @@ for the JavaScript code in this tag.
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org15c3fe1">1. Setup</a></li>
<li><a href="#org10e7b57">2. Signal Processing</a>
<li><a href="#orgeb3ab60">1. Experimental Setup</a></li>
<li><a href="#org17b1fa2">2. Signal Processing</a>
<ul>
<li><a href="#orgce9374f">2.1. Load data</a></li>
<li><a href="#orgf5bde3c">2.2. Time Domain Data</a></li>
<li><a href="#org53fa3f2">2.3. Computation of the ASD of the measured voltage</a></li>
<li><a href="#orgdb1374d">2.4. Scaling to take into account the sensibility of the geophone and the voltage amplifier</a></li>
<li><a href="#org07e8527">2.5. Computation of the ASD of the velocity</a></li>
<li><a href="#orgf3d4fe6">2.6. Transfer function between the two geophones</a></li>
<li><a href="#org1336e7b">2.7. Estimation of the sensor noise</a></li>
<li><a href="#orge77a4ee">2.1. Load data</a></li>
<li><a href="#org6b7befd">2.2. Time Domain Data</a></li>
<li><a href="#org4e41681">2.3. Computation of the ASD of the measured voltage</a></li>
<li><a href="#org83cd67f">2.4. Scaling to take into account the sensibility of the geophone and the voltage amplifier</a></li>
<li><a href="#orge0f93a9">2.5. Computation of the ASD of the velocity</a></li>
<li><a href="#orgd178821">2.6. Transfer function between the two geophones</a></li>
<li><a href="#orgaf0dac1">2.7. Estimation of the sensor noise</a></li>
</ul>
</li>
<li><a href="#orgaee4ca7">3. Compare axis</a>
<li><a href="#org126e178">3. Compare axis</a>
<ul>
<li><a href="#org4a91ee9">3.1. Load data</a></li>
<li><a href="#org36adfc5">3.2. Compare PSD</a></li>
<li><a href="#org698e322">3.3. Compare TF</a></li>
<li><a href="#org3f6cac8">3.1. Load data</a></li>
<li><a href="#org608de0c">3.2. Compare PSD</a></li>
<li><a href="#org4db3872">3.3. Compare TF</a></li>
</ul>
</li>
<li><a href="#orge3ff149">4. Appendix</a>
<li><a href="#org66ca70f">4. Appendix</a>
<ul>
<li><a href="#org2060421">4.1. Computation of coherence from PSD and CSD</a></li>
<li><a href="#org7492edb">4.1. Computation of coherence from PSD and CSD</a></li>
</ul>
</li>
</ul>
</div>
</div>
<div id="outline-container-org15c3fe1" class="outline-2">
<h2 id="org15c3fe1"><span class="section-number-2">1</span> Setup</h2>
<div id="outline-container-orgeb3ab60" class="outline-2">
<h2 id="orgeb3ab60"><span class="section-number-2">1</span> Experimental Setup</h2>
<div class="outline-text-2" id="text-1">
<p>
Two L22 geophones are used.
@ -319,14 +319,14 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH
</p>
<div id="org2286cc3" class="figure">
<div id="orgfe0dca3" class="figure">
<p><img src="./figs/setup.jpg" alt="setup.jpg" width="500px" />
</p>
<p><span class="figure-number">Figure 1: </span>Setup</p>
</div>
<div id="org9446018" class="figure">
<div id="orgbcec6a9" class="figure">
<p><img src="./figs/geophones.jpg" alt="geophones.jpg" width="500px" />
</p>
<p><span class="figure-number">Figure 2: </span>Geophones</p>
@ -334,16 +334,22 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH
</div>
</div>
<div id="outline-container-org10e7b57" class="outline-2">
<h2 id="org10e7b57"><span class="section-number-2">2</span> Signal Processing</h2>
<div id="outline-container-org17b1fa2" class="outline-2">
<h2 id="org17b1fa2"><span class="section-number-2">2</span> Signal Processing</h2>
<div class="outline-text-2" id="text-2">
<p>
The Matlab computing file for this part is accessible <a href="signal_processing.m">here</a>.
The <code>mat</code> file containing the measurement data is accessible <a href="mat/data_001.mat">here</a>.
</p>
</div>
<div id="outline-container-orgce9374f" class="outline-3">
<h3 id="orgce9374f"><span class="section-number-3">2.1</span> Load data</h3>
<div id="outline-container-orge77a4ee" class="outline-3">
<h3 id="orge77a4ee"><span class="section-number-3">2.1</span> Load data</h3>
<div class="outline-text-3" id="text-2-1">
<p>
We load the data of the z axis of two geophones.
</p>
<div class="org-src-container">
<pre class="src src-matlab">load<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-string">'mat/data_001.mat', 't', 'x1', 'x2'</span><span class="org-rainbow-delimiters-depth-1">)</span>;
dt = t<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-highlight-numbers-number">2</span><span class="org-rainbow-delimiters-depth-1">)</span> <span class="org-type">-</span> t<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-highlight-numbers-number">1</span><span class="org-rainbow-delimiters-depth-1">)</span>;
@ -352,8 +358,8 @@ dt = t<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-high
</div>
</div>
<div id="outline-container-orgf5bde3c" class="outline-3">
<h3 id="orgf5bde3c"><span class="section-number-3">2.2</span> Time Domain Data</h3>
<div id="outline-container-org6b7befd" class="outline-3">
<h3 id="org6b7befd"><span class="section-number-3">2.2</span> Time Domain Data</h3>
<div class="outline-text-3" id="text-2-2">
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-type">figure</span>;
@ -368,7 +374,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="org2d30e99" class="figure">
<div id="org0161433" class="figure">
<p><img src="figs/data_time_domain.png" alt="data_time_domain.png" />
</p>
<p><span class="figure-number">Figure 3: </span>Time domain Data</p>
@ -388,7 +394,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="org8e28356" class="figure">
<div id="orgab64b52" class="figure">
<p><img src="figs/data_time_domain_zoom.png" alt="data_time_domain_zoom.png" />
</p>
<p><span class="figure-number">Figure 4: </span>Time domain Data - Zoom</p>
@ -396,8 +402,8 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
</div>
<div id="outline-container-org53fa3f2" class="outline-3">
<h3 id="org53fa3f2"><span class="section-number-3">2.3</span> Computation of the ASD of the measured voltage</h3>
<div id="outline-container-org4e41681" class="outline-3">
<h3 id="org4e41681"><span class="section-number-3">2.3</span> Computation of the ASD of the measured voltage</h3>
<div class="outline-text-3" id="text-2-3">
<p>
We first define the parameters for the frequency domain analysis.
@ -416,12 +422,12 @@ Fs = <span class="org-highlight-numbers-number">1</span><span class="org-type">/
</div>
</div>
<div id="outline-container-orgdb1374d" class="outline-3">
<h3 id="orgdb1374d"><span class="section-number-3">2.4</span> Scaling to take into account the sensibility of the geophone and the voltage amplifier</h3>
<div id="outline-container-org83cd67f" class="outline-3">
<h3 id="org83cd67f"><span class="section-number-3">2.4</span> Scaling to take into account the sensibility of the geophone and the voltage amplifier</h3>
<div class="outline-text-3" id="text-2-4">
<p>
The Geophone used are L22.
Their sensibility are shown on figure <a href="#org59aaa6d">5</a>.
Their sensibility are shown on figure <a href="#orgf03c4b4">5</a>.
</p>
<div class="org-src-container">
@ -432,7 +438,7 @@ S = <span class="org-rainbow-delimiters-depth-1">(</span>s<span class="org-type"
</div>
<div id="org59aaa6d" class="figure">
<div id="orgf03c4b4" class="figure">
<p><img src="figs/geophone_sensibility.png" alt="geophone_sensibility.png" />
</p>
<p><span class="figure-number">Figure 5: </span>Sensibility of the Geophone</p>
@ -463,11 +469,11 @@ We further divide the result by the sensibility of the Geophone to obtain the AS
</div>
</div>
<div id="outline-container-org07e8527" class="outline-3">
<h3 id="org07e8527"><span class="section-number-3">2.5</span> Computation of the ASD of the velocity</h3>
<div id="outline-container-orge0f93a9" class="outline-3">
<h3 id="orge0f93a9"><span class="section-number-3">2.5</span> Computation of the ASD of the velocity</h3>
<div class="outline-text-3" id="text-2-5">
<p>
The ASD of the measured velocity is shown on figure <a href="#orgd1800a0">6</a>.
The ASD of the measured velocity is shown on figure <a href="#org9a56511">6</a>.
</p>
<div class="org-src-container">
@ -484,7 +490,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="orgd1800a0" class="figure">
<div id="org9a56511" class="figure">
<p><img src="figs/psd_velocity.png" alt="psd_velocity.png" />
</p>
<p><span class="figure-number">Figure 6: </span>Spectral density of the velocity</p>
@ -492,16 +498,16 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
</div>
<div id="outline-container-orgf3d4fe6" class="outline-3">
<h3 id="orgf3d4fe6"><span class="section-number-3">2.6</span> Transfer function between the two geophones</h3>
<div id="outline-container-orgd178821" class="outline-3">
<h3 id="orgd178821"><span class="section-number-3">2.6</span> Transfer function between the two geophones</h3>
<div class="outline-text-3" id="text-2-6">
<p>
We here compute the transfer function from one geophone to the other.
The result is shown on figure <a href="#org6abd5e1">7</a>.
The result is shown on figure <a href="#orgb6c07f9">7</a>.
</p>
<p>
We also compute the coherence between the two signals (figure <a href="#org80d61f1">8</a>).
We also compute the coherence between the two signals (figure <a href="#org1b45a36">8</a>).
</p>
<div class="org-src-container">
@ -510,7 +516,7 @@ We also compute the coherence between the two signals (figure <a href="#org80d61
</div>
<div id="org6abd5e1" class="figure">
<div id="orgb6c07f9" class="figure">
<p><img src="figs/tf_geophones.png" alt="tf_geophones.png" />
</p>
<p><span class="figure-number">Figure 7: </span>Estimated transfer function between the two geophones</p>
@ -522,7 +528,7 @@ We also compute the coherence between the two signals (figure <a href="#org80d61
</div>
<div id="org80d61f1" class="figure">
<div id="org1b45a36" class="figure">
<p><img src="figs/coh_geophones.png" alt="coh_geophones.png" />
</p>
<p><span class="figure-number">Figure 8: </span>Cohererence between the signals of the two geophones</p>
@ -530,8 +536,8 @@ We also compute the coherence between the two signals (figure <a href="#org80d61
</div>
</div>
<div id="outline-container-org1336e7b" class="outline-3">
<h3 id="org1336e7b"><span class="section-number-3">2.7</span> Estimation of the sensor noise</h3>
<div id="outline-container-orgaf0dac1" class="outline-3">
<h3 id="orgaf0dac1"><span class="section-number-3">2.7</span> Estimation of the sensor noise</h3>
<div class="outline-text-3" id="text-2-7">
<p>
The technique to estimate the sensor noise is taken from <a class='org-ref-reference' href="#barzilai98_techn_measur_noise_sensor_presen">barzilai98_techn_measur_noise_sensor_presen</a>.
@ -561,11 +567,11 @@ where:
</ul>
<p>
The <code>mscohere</code> function is compared with this formula on Appendix (section <a href="#org72623b2">4.1</a>), it is shown that it is identical.
The <code>mscohere</code> function is compared with this formula on Appendix (section <a href="#orgd085438">4.1</a>), it is shown that it is identical.
</p>
<p>
Figure <a href="#org88b4597">9</a> illustrate a block diagram model of the system used to determine the sensor noise of the geophone.
Figure <a href="#orga4af110">9</a> illustrate a block diagram model of the system used to determine the sensor noise of the geophone.
</p>
<p>
@ -577,7 +583,7 @@ Each sensor has noise \(N\) and \(M\).
</p>
<div id="org88b4597" class="figure">
<div id="orga4af110" class="figure">
<p><img src="figs/huddle-test.png" alt="huddle-test.png" />
</p>
<p><span class="figure-number">Figure 9: </span>Huddle test block diagram</p>
@ -592,7 +598,7 @@ We also assume that \(H_1 = H_2 = 1\).
We then obtain:
</p>
\begin{equation}
\label{org5b4a541}
\label{org65b3ddf}
\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}
@ -600,23 +606,23 @@ We then obtain:
Since the input signal \(U\) and the instrumental noise \(N\) are incoherent:
</p>
\begin{equation}
\label{orga04bdf3}
\label{org14038cd}
|G_X(\omega)| = |G_N(\omega)| + |G_U(\omega)|
\end{equation}
<p>
From equations \eqref{org5b4a541} and \eqref{orga04bdf3}, we finally obtain
From equations \eqref{org65b3ddf} and \eqref{org14038cd}, we finally obtain
</p>
<div class="important">
\begin{equation}
\label{orga57b1ae}
\label{org84c455c}
|G_N(\omega)| = |G_X(\omega)| \left( 1 - \sqrt{\gamma_{XY}^2(\omega)} \right)
\end{equation}
</div>
<p>
The instrumental noise is computed below. The result in V<sup>2</sup>/Hz is shown on figure <a href="#orgd3329eb">10</a>.
The instrumental noise is computed below. The result in V<sup>2</sup>/Hz is shown on figure <a href="#orgd0b903d">10</a>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">pxxN = pxx1<span class="org-type">.*</span><span class="org-rainbow-delimiters-depth-1">(</span><span class="org-highlight-numbers-number">1</span> <span class="org-type">-</span> coh12<span class="org-rainbow-delimiters-depth-1">)</span>;
@ -637,14 +643,14 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="orgd3329eb" class="figure">
<div id="orgd0b903d" class="figure">
<p><img src="figs/intrumental_noise_V.png" alt="intrumental_noise_V.png" />
</p>
<p><span class="figure-number">Figure 10: </span>Instrumental Noise and Measurement in \(V^2/Hz\)</p>
</div>
<p>
This is then further converted into velocity and compared with the ground velocity measurement. (figure <a href="#org6fa2e55">11</a>)
This is then further converted into velocity and compared with the ground velocity measurement. (figure <a href="#org0542ab9">11</a>)
</p>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-type">figure</span>;
@ -660,7 +666,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="org6fa2e55" class="figure">
<div id="org0542ab9" class="figure">
<p><img src="figs/intrumental_noise_velocity.png" alt="intrumental_noise_velocity.png" />
</p>
<p><span class="figure-number">Figure 11: </span>Instrumental Noise and Measurement in \(m/s/\sqrt{Hz}\)</p>
@ -669,13 +675,26 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
</div>
<div id="outline-container-orgaee4ca7" class="outline-2">
<h2 id="orgaee4ca7"><span class="section-number-2">3</span> Compare axis</h2>
<div id="outline-container-org126e178" class="outline-2">
<h2 id="org126e178"><span class="section-number-2">3</span> Compare axis</h2>
<div class="outline-text-2" id="text-3">
<p>
The Matlab computing file for this part is accessible <a href="compare_axis.m">here</a>.
The <code>mat</code> files containing the measurement data are accessible with the following links:
</p>
<ul class="org-ul">
<li>z axis: <a href="mat/data_001.mat">here</a>.</li>
<li>east axis: <a href="mat/data_002.mat">here</a>.</li>
<li>north axis: <a href="mat/data_003.mat">here</a>.</li>
</ul>
</div>
<div id="outline-container-org4a91ee9" class="outline-3">
<h3 id="org4a91ee9"><span class="section-number-3">3.1</span> Load data</h3>
<div id="outline-container-org3f6cac8" class="outline-3">
<h3 id="org3f6cac8"><span class="section-number-3">3.1</span> Load data</h3>
<div class="outline-text-3" id="text-3-1">
<p>
We first load the data for the three axis.
</p>
<div class="org-src-container">
<pre class="src src-matlab">z = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-string">'mat/data_001.mat', 't', 'x1', 'x2'</span><span class="org-rainbow-delimiters-depth-1">)</span>;
east = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-string">'mat/data_002.mat', 't', 'x1', 'x2'</span><span class="org-rainbow-delimiters-depth-1">)</span>;
@ -685,9 +704,12 @@ north = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="or
</div>
</div>
<div id="outline-container-org36adfc5" class="outline-3">
<h3 id="org36adfc5"><span class="section-number-3">3.2</span> Compare PSD</h3>
<div id="outline-container-org608de0c" class="outline-3">
<h3 id="org608de0c"><span class="section-number-3">3.2</span> Compare PSD</h3>
<div class="outline-text-3" id="text-3-2">
<p>
The PSD for each axis of the two geophones are computed.
</p>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-rainbow-delimiters-depth-1">[</span>pz1, fz<span class="org-rainbow-delimiters-depth-1">]</span> = pwelch<span class="org-rainbow-delimiters-depth-1">(</span>z.x1, hanning<span class="org-rainbow-delimiters-depth-2">(</span>ceil<span class="org-rainbow-delimiters-depth-3">(</span>length<span class="org-rainbow-delimiters-depth-4">(</span>z.x1<span class="org-rainbow-delimiters-depth-4">)</span><span class="org-type">/</span><span class="org-highlight-numbers-number">100</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-highlight-numbers-number">1</span><span class="org-type">/</span><span class="org-rainbow-delimiters-depth-2">(</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">2</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-type">-</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">1</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span><span class="org-rainbow-delimiters-depth-1">)</span>;
<span class="org-rainbow-delimiters-depth-1">[</span>pz2, <span class="org-type">~</span><span class="org-rainbow-delimiters-depth-1">]</span> = pwelch<span class="org-rainbow-delimiters-depth-1">(</span>z.x2, hanning<span class="org-rainbow-delimiters-depth-2">(</span>ceil<span class="org-rainbow-delimiters-depth-3">(</span>length<span class="org-rainbow-delimiters-depth-4">(</span>z.x2<span class="org-rainbow-delimiters-depth-4">)</span><span class="org-type">/</span><span class="org-highlight-numbers-number">100</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-highlight-numbers-number">1</span><span class="org-type">/</span><span class="org-rainbow-delimiters-depth-2">(</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">2</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-type">-</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">1</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span><span class="org-rainbow-delimiters-depth-1">)</span>;
@ -700,8 +722,11 @@ north = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="or
</pre>
</div>
<p>
We compare them. The result is shown on figure <a href="#orge0ebe78">12</a>.
</p>
<div id="orgc2c3d69" class="figure">
<div id="orge0ebe78" class="figure">
<p><img src="figs/compare_axis_psd.png" alt="compare_axis_psd.png" />
</p>
<p><span class="figure-number">Figure 12: </span>Compare the measure PSD of the two geophones for the three axis</p>
@ -709,9 +734,14 @@ north = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="or
</div>
</div>
<div id="outline-container-org698e322" class="outline-3">
<h3 id="org698e322"><span class="section-number-3">3.3</span> Compare TF</h3>
<div id="outline-container-org4db3872" class="outline-3">
<h3 id="org4db3872"><span class="section-number-3">3.3</span> Compare TF</h3>
<div class="outline-text-3" id="text-3-3">
<p>
The transfer functions from one geophone to the other are also computed for each axis.
The result is shown on figure <a href="#org2a4c622">13</a>.
</p>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-rainbow-delimiters-depth-1">[</span>Tz, fz<span class="org-rainbow-delimiters-depth-1">]</span> = tfestimate<span class="org-rainbow-delimiters-depth-1">(</span>z.x1, z.x2, hanning<span class="org-rainbow-delimiters-depth-2">(</span>ceil<span class="org-rainbow-delimiters-depth-3">(</span>length<span class="org-rainbow-delimiters-depth-4">(</span>z.x1<span class="org-rainbow-delimiters-depth-4">)</span><span class="org-type">/</span><span class="org-highlight-numbers-number">100</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-highlight-numbers-number">1</span><span class="org-type">/</span><span class="org-rainbow-delimiters-depth-2">(</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">2</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-type">-</span>z.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">1</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span><span class="org-rainbow-delimiters-depth-1">)</span>;
<span class="org-rainbow-delimiters-depth-1">[</span>Te, fe<span class="org-rainbow-delimiters-depth-1">]</span> = tfestimate<span class="org-rainbow-delimiters-depth-1">(</span>east.x1, east.x2, hanning<span class="org-rainbow-delimiters-depth-2">(</span>ceil<span class="org-rainbow-delimiters-depth-3">(</span>length<span class="org-rainbow-delimiters-depth-4">(</span>east.x1<span class="org-rainbow-delimiters-depth-4">)</span><span class="org-type">/</span><span class="org-highlight-numbers-number">100</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-rainbow-delimiters-depth-2">[]</span>, <span class="org-highlight-numbers-number">1</span><span class="org-type">/</span><span class="org-rainbow-delimiters-depth-2">(</span>east.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">2</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-type">-</span>east.t<span class="org-rainbow-delimiters-depth-3">(</span><span class="org-highlight-numbers-number">1</span><span class="org-rainbow-delimiters-depth-3">)</span><span class="org-rainbow-delimiters-depth-2">)</span><span class="org-rainbow-delimiters-depth-1">)</span>;
@ -720,7 +750,7 @@ north = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="or
</div>
<div id="orgf7b4d80" class="figure">
<div id="org2a4c622" class="figure">
<p><img src="figs/compare_tf_axis.png" alt="compare_tf_axis.png" />
</p>
<p><span class="figure-number">Figure 13: </span>Compare the transfer function from one geophone to the other for the 3 axis</p>
@ -728,15 +758,16 @@ north = load<span class="org-rainbow-delimiters-depth-1">(</span><span class="or
</div>
</div>
</div>
<div id="outline-container-orge3ff149" class="outline-2">
<h2 id="orge3ff149"><span class="section-number-2">4</span> Appendix</h2>
<div id="outline-container-org66ca70f" class="outline-2">
<h2 id="org66ca70f"><span class="section-number-2">4</span> Appendix</h2>
<div class="outline-text-2" id="text-4">
</div>
<div id="outline-container-org2060421" class="outline-3">
<h3 id="org2060421"><span class="section-number-3">4.1</span> Computation of coherence from PSD and CSD</h3>
<div id="outline-container-org7492edb" class="outline-3">
<h3 id="org7492edb"><span class="section-number-3">4.1</span> Computation of coherence from PSD and CSD</h3>
<div class="outline-text-3" id="text-4-1">
<p>
<a id="org72623b2"></a>
<a id="orgd085438"></a>
</p>
<div class="org-src-container">
<pre class="src src-matlab">load<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-string">'mat/data_001.mat', 't', 'x1', 'x2'</span><span class="org-rainbow-delimiters-depth-1">)</span>;
@ -767,7 +798,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="org3184d2c" class="figure">
<div id="org0bcdbf7" class="figure">
<p><img src="figs/comp_coherence_formula.png" alt="comp_coherence_formula.png" />
</p>
<p><span class="figure-number">Figure 14: </span>Comparison of <code>mscohere</code> and manual computation</p>
@ -785,7 +816,7 @@ xlim<span class="org-rainbow-delimiters-depth-1">(</span><span class="org-rainbo
</div>
<div id="postamble" class="status">
<p class="author">Author: Thomas Dehaeze</p>
<p class="date">Created: 2019-04-18 jeu. 17:02</p>
<p class="date">Created: 2019-04-18 jeu. 17:11</p>
<p class="validation"><a href="http://validator.w3.org/check?uri=referer">Validate</a></p>
</div>
</body>

29
huddle-test-geophones/index.org
View File

@ -18,7 +18,7 @@
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:
* Setup
* Experimental Setup
Two L22 geophones are used.
They are placed on the ID31 granite.
They are leveled.
@ -37,6 +37,13 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH
[[file:./figs/geophones.jpg]]
* Signal Processing
:PROPERTIES:
:header-args:matlab+: :tangle signal_processing.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
The Matlab computing file for this part is accessible [[file:signal_processing.m][here]].
The =mat= file containing the measurement data is accessible [[file:mat/data_001.mat][here]].
** Matlab Init :noexport:ignore:
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
@ -44,6 +51,7 @@ The voltage amplifiers include a low pass filter with a cut-off frequency at 1kH
** 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);
@ -346,12 +354,23 @@ This is then further converted into velocity and compared with the ground veloci
[[file:figs/intrumental_noise_velocity.png]]
* Compare axis
:PROPERTIES:
:header-args:matlab+: :tangle compare_axis.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
The Matlab computing file for this part is accessible [[file:compare_axis.m][here]].
The =mat= files containing the measurement data are accessible with the following links:
- z axis: [[file:mat/data_001.mat][here]].
- east axis: [[file:mat/data_002.mat][here]].
- north axis: [[file:mat/data_003.mat][here]].
** Matlab Init :noexport:ignore:
#+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');
@ -359,6 +378,7 @@ This is then further converted into velocity and compared with the ground veloci
#+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)));
@ -370,7 +390,8 @@ This is then further converted into velocity and compared with the ground veloci
[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
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');
@ -398,6 +419,9 @@ This is then further converted into velocity and compared with the ground veloci
[[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)));
@ -442,6 +466,7 @@ This is then further converted into velocity and compared with the ground veloci
#+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>>

233
huddle-test-geophones/signal_processing.m Normal file
View File

@ -0,0 +1,233 @@
% Matlab Init :noexport:ignore:
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
%% Initialize ans with org-babel
ans = 0;
% Load data
% We load the data of the z axis of two geophones.
load('mat/data_001.mat', 't', 'x1', 'x2');
dt = t(2) - t(1);
% Time Domain Data
figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([t(1), t(end)]);
% #+NAME: fig:data_time_domain
% #+CAPTION: Time domain Data
% #+RESULTS: fig:data_time_domain
% [[file:figs/data_time_domain.png]]
figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([0 1]);
% Computation of the ASD of the measured voltage
% We first define the parameters for the frequency domain analysis.
win = hanning(ceil(length(x1)/100));
Fs = 1/dt;
[pxx1, f] = pwelch(x1, win, [], [], Fs);
[pxx2, ~] = pwelch(x2, win, [], [], Fs);
% 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]].
S0 = 88; % Sensitivity [V/(m/s)]
f0 = 2; % Cut-off frequnecy [Hz]
S = (s/2/pi/f0)/(1+s/2/pi/f0);
figure;
bodeFig({S});
ylabel('Amplitude [V/(m/s)]')
% #+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.
G0 = 60; % [dB]
G = G0/(1+s/2/pi/1000);
% 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}$.
scaling = 1./squeeze(abs(freqresp(G, f, 'Hz')))./squeeze(abs(freqresp(S, f, 'Hz')));
% Computation of the ASD of the velocity
% The ASD of the measured velocity is shown on figure [[fig:psd_velocity]].
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]);
% 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]]).
[T12, ~] = tfestimate(x1, x2, win, [], [], Fs);
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]);
% #+NAME: fig:tf_geophones
% #+CAPTION: Estimated transfer function between the two geophones
% #+RESULTS: fig:tf_geophones
% [[file:figs/tf_geophones.png]]
[coh12, ~] = mscohere(x1, x2, win, [], [], Fs);
figure;
plot(f, coh12);
set(gca, 'xscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Coherence');
ylim([0,1]); xlim([1, 500]);
% 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]].
pxxN = pxx1.*(1 - coh12);
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]);
% #+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]])
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]);