This commit is contained in:
Thomas Dehaeze 2019-05-10 09:13:09 +02:00
commit 16c222c229
184 changed files with 28668 additions and 20 deletions

46
.gitignore vendored
View File

@ -2,3 +2,49 @@ auto/
*.tex
**/figs/*.pdf
**/figs/*.svg
=======
# Emacs
auto/
# Simulink Real Time
*bio.m
*pt.m
*ref.m
*ri.m
*xcp.m
*.mldatx
*.slxc
*.xml
*_slrt_rtw/
# data
data/
# Windows default autosave extension
*.asv
# OSX / *nix default autosave extension
*.m~
# Compiled MEX binaries (all platforms)
*.mex*
# Packaged app and toolbox files
*.mlappinstall
*.mltbx
# Generated helpsearch folders
helpsearch*/
# Simulink code generation folders
slprj/
sccprj/
# Matlab code generation folders
codegen/
# Simulink autosave extension
*.autosave
# Octave session info
octave-workspace

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1,196 @@
#+TITLE:Measurement of the sample vibrations when rotating the Spindle
: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:
* Experimental Setup
* Signal Processing
** Matlab Init :noexport:ignore:
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
** Load Data
Measurement =data_001.mat= corresponds to a measurement where the spindle is not turning and =data_002.mat= where the spindle is turning at 1rpm.
=x1= is the signal coming from the geophone located on the marble and =x2= is the signal from the geophone located on the sample station.
#+begin_src matlab :results none
data1 = load('mat/data_001.mat', 't', 'x1', 'x2');
data2 = load('mat/data_002.mat', 't', 'x1', 'x2');
#+end_src
** Time domain Data
#+begin_src matlab :results none
figure;
hold on;
plot(data1.t, data1.x1);
plot(data2.t, data2.x1);
hold off;
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(data1.t, data1.x2);
plot(data2.t, data2.x2)
hold off;
#+end_src
** ASD and Frequency domain data
#+begin_src matlab :results none
dt = data1.t(2) - data1.t(1);
Fs = 1/dt;
windows_psd = hanning(ceil(10*Fs));
#+end_src
#+begin_src matlab :results none
[pxx1m, f] = pwelch(data1.x1, windows_psd, [], [], Fs); f(1) = []; pxx1m(1) = [];
[pxx1h, ~] = pwelch(data1.x2, windows_psd, [], [], Fs); pxx1h(1) = [];
[pxx2m, ~] = pwelch(data2.x1, windows_psd, [], [], Fs); pxx2m(1) = [];
[pxx2h, ~] = pwelch(data2.x2, windows_psd, [], [], Fs); pxx2h(1) = [];
#+end_src
** Some plots
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(pxx1m));
plot(f, sqrt(pxx2m));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(pxx1h));
plot(f, sqrt(pxx2h));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(pxx2m));
plot(f, sqrt(pxx2h));
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, cumtrapz(f, pxx1m))
plot(f, cumtrapz(f, pxx2m))
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('CAS [m]')
#+end_src
** 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
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^(60/G0_db); % [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'))); scaling(1) = 0;
#+end_src
** Computation of the ASD of the velocity
#+begin_src matlab :results none
figure;
hold on;
plot(f, sqrt(pxx1h).*scaling);
plot(f, sqrt(pxx2h).*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
#+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
** RMS value of the difference between the two geophones
We also compute the Power Spectral Density of the difference between the two geophones. This is done in order to estimate the relative displacement of the sample with respect to the granite.
#+begin_src matlab :results none
[pxxd1, ~] = pwelch(data1.x2-data1.x1, windows_psd, [], [], Fs); pxxd1(1) = [];
[pxxd2, ~] = pwelch(data2.x2-data2.x1, windows_psd, [], [], Fs); pxxd2(1) = [];
#+end_src
#+begin_src matlab :results none
figure;
hold on;
plot(f, (sqrt(pxxd1).*scaling)./(2*pi*f));
plot(f, (sqrt(pxxd2).*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
#+begin_src matlab :results none
psd_d1 = ((sqrt(pxxd1).*scaling)./(2*pi*f)).^2;
psd_d2 = ((sqrt(pxxd2).*scaling)./(2*pi*f)).^2;
df = f(2) - f(1);
figure;
hold on;
plot(f, sqrt(cumsum(df.*psd_d1, 'reverse')));
plot(f, sqrt(cumsum(df.*psd_d2, 'reverse')));
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('CAS $\left[m\right]$')
xlim([0.1, 500]);
#+end_src

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1,12 @@
#+TITLE:Measurement of the sample vibrations when rotating the Spindle
Measurement:
- one geophone on the marble
- one geophone at the sample position
Every stage is powered on.
| Data file | Description |
|--------------+----------------------------|
| data_001.mat | nothing is turning |
| data_002.mat | spindle is turning at 1rpm |

View File

@ -0,0 +1,53 @@
%%
Tsim = 100; % [s]
%%
tg = slrt;
%% TODO - Build this application if updated
%%
if tg.Connected == "Yes"
if tg.Status == "running"
disp('Target is Running, Stopping...');
tg.stop;
while tg.Status == "running"
pause(1);
end
disp('Target is Stopped');
end
if tg.Status == "stopped"
disp('Load the Application');
tg.load('disturbance_measurement');
%% Run the application
disp('Starting the Application');
tg.start;
pause(Tsim);
tg.stop;
end
else
error("The target computer is not connected");
end
%%
f = SimulinkRealTime.openFTP(tg);
cd(f, 'data/disturbance_measurement/');
mget(f, 'data_001.dat', 'data');
close(f);
data = SimulinkRealTime.utils.getFileScopeData('data/data_001.dat').data;
%%
n = 17;
while isfile(['mat/data_', num2str(n, '%03d'), '.mat'])
disp('File exists.');
if input(['Are you sure you want to override the file ', 'mat/data_', ...
num2str(n, '%03d'), '.mat', ' ? [Y/n]']) == 'Y'
break;
end
n = input('What should be the measurement number?');
end
save(['mat/data_', num2str(n, '%03d'), '.mat'], 'data');

View File

Binary file not shown.

BIN
equipment/equipment.html Normal file

Binary file not shown.

28
equipment/equipment.org Normal file
View File

@ -0,0 +1,28 @@
#+TITLE: Equipment used to make the measurements
* Geophone
L22
* Speedgoat
ADC and DAC:
* Voltage amplifier
#+name: fig:voltage_amplifier
#+caption: Picture of the voltage amplifier DLPVA
#+attr_latex: :scale 1
[[file:./img/DLPVA_W_R2.jpg]]
*Specifications*:
- Switchable gain up to 100 dB (x 100,000)
- Bandwidth DC to 100 kHz
- DC-drift 0.6 µV/°C
- Input noise down to 0.4 nV/√Hz
- Switchable AC/DC-coupling
- Local and remote control
- Input impedance up to 1 TΩ
The documentation of the voltage amplifier is accessible [[file:documents/dlpva_datasheet.pdf][here]].
The exact model is DLPVA-100-B-D

Binary file not shown.

After

Width:  |  Height:  |  Size: 36 KiB

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]);

3
huddle-test-geophones/figs/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
*.svg
*.pdf
*.tex

Binary file not shown.

After

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 134 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 65 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 180 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 42 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 47 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.7 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.3 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 169 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 153 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 136 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.8 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 99 KiB

Binary file not shown.

Binary file not shown.

View File

@ -0,0 +1,561 @@
#+TITLE:Huddle Test of the L22 Geophones
: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:
* 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 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
: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>>
#+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^(60/G0_db); % [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 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');
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

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -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},
}

View File

@ -0,0 +1,47 @@
tg = slrt;
%% TODO - Build this application if updated
%%
if tg.Connected == "Yes"
if tg.Status == "stopped"
%% Load the application
tg.load('test');
%% Run the application
tg.start;
pause(10);
tg.stop;
%% Load the data
f = SimulinkRealTime.openFTP(tg);
mget(f, 'data/data_001.dat');
close(f);
end
end
%% Convert the Data
data = SimulinkRealTime.utils.getFileScopeData('data/data_001.dat').data;
t = data(:, end);
x1 = data(:, 1);
x2 = data(:, 2);
save('mat/data_003.mat', 't', 'x1', 'x2');
%% Plot the data
figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off
xlabel('Time [s]');
ylabel('Voltage [V]');
%% Compute the PSD
dt = t(2)-t(1);
window_L = ceil(length(x1)/10);
window_han = .5*(1 - cos(2*pi*(1:window_L)'/(window_L+1)));
[pxx, f] = pwelch(x1, window_han, 0, [], 1/dt);

View File

@ -0,0 +1 @@
Ts = 1e-3; % [s]

View File

@ -0,0 +1,252 @@
% 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 = 10^(G0/20)/(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*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]);
% #+NAME: fig:psd_velocity
% #+CAPTION: 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]]);
figure;
hold on;
plot(f, (pxx1.*scaling./f).^2);
plot(f, (pxx2.*scaling./f).^2);
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]);

Binary file not shown.

View File

@ -2,26 +2,13 @@
:drawer:
#+STARTUP: overview
#+OPTIONS: toc:2
#+OPTIONS: num:2
#+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: <script src="js/jquery.min.js"></script>
#+HTML_HEAD: <script src="js/bootstrap.min.js"></script>
#+HTML_HEAD: <script src="js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script src="js/readtheorg.js"></script>
#+LATEX_CLASS: cleanreport
#+LaTeX_CLASS_OPTIONS: [tocnp, secbreak, minted]
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs
#+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>
:end:
* Measurements of the dynamics of the station
@ -141,3 +128,10 @@ The goal is to estimate all the error motions induced by the Spindle
* Ressources
[[file:actuators-sensors/index.org][Actuators and Sensors]]
* Other measurements
- [[file:huddle-test-geophones/index.org][Huddle Test - Geophones]]
- [[file:disturbance-measurement/index.org][Disturbance Measurement]]
- [[file:slip-ring-test/index.org][Slip Ring - Noise measurement]]
- [[file:static-measurements/index.org][Control System Measurement]]
- [[file:equipment/equipment.org][Equipment used for the measurements]]

86
readme.org Normal file
View File

@ -0,0 +1,86 @@
* Usefull Simulink Real Time commands
* Get data from the target computer
We first copy the =dat= file from the target computer to the host computer:
#+begin_src matlab
tg=slrt;
f=SimulinkRealTime.openFTP(tg);
mget(f, 'DATA_001.dat', 'local_folder');
close(f);
#+end_src
We then import the =dat= file to the workspace:
#+begin_src matlab
matlab_data = SimulinkRealTime.utils.getFileScopeData('DATA_001.dat');
#+end_src
* Commands with the target object
https://fr.mathworks.com/help/xpc/api/simulinkrealtime.target.html
| Connect | No Yes |
| Status | stopped runing |
| start | Start execution of real-time application on target computer |
| stop | Stop execution of real-time application on target computer |
| tg.viewTargetScreen | Show target computer screen |
| ping | Test communication between development and target computers |
| reboot | Restart target computer |
| close | Close connection between development and target computers |
| load | Download real-time application to target computer |
| unload | Remove real-time application from target computer |
| addscope | Create a scope of specified type |
| getscope | Return scope identified by scope number |
| remscope | Remove scope from target computer |
| getlog | Portion of output logs from target object |
| importLogData | Import buffered logging data to the active session of the Simulation Data Inspector |
| getsignal | Value of signal |
| getsignalid | Signal index from signal hierarchical name |
| getsignalidsfromlabel | Vector of signal indices |
| getsignallabel | Signal label for signal index |
| getsignalname | Signal name from index list |
| getparam | Read value of observable parameter in real-time application |
| setparam | Change value of tunable parameter in real-time application |
| getparamid | Parameter index from parameter hierarchical name |
| getparamname | Block path and parameter name from parameter index |
| loadparamset | Restore parameter values saved in specified file |
| saveparamset | Save real-time application parameter values |
| startProfiler | Start profiling service on target computer |
| stopProfiler | Stop profiling service on target computer |
| getProfilerData | Retrieve profile data object |
| resetProfiler | Reset profiling service state to Ready |
| getDiskSpace | Return free space and total space on the drive, in bytes |
* FTP access to the target computer
https://fr.mathworks.com/help/xpc/api/simulinkrealtime.openftp.html?s_tid=doc_ta
First run the following commands to have the =FTP= Object:
#+begin_src matlab
tg=slrt;
f=SimulinkRealTime.openFTP(tg);
#+end_src
Then, the =f= object can be used to access the filesystem on the target computer.
| cd | | |
| dir | | |
| mget | Used to download data from the target host | =f.mget('data.dat', 'local_folder')= |
| mkdir | | |
| mput | | |
| rename | | |
| rmdir | | |
| close | | |
* ELMO
tutorials: https://www.elmomc.com/products/application-studio/easii/easii-tutorials/
* Low Pass Filter
R = 1KOhm
C = 1muF
Fc = 1kHz

11
slip-ring-test/analysis.m Normal file
View File

@ -0,0 +1,11 @@
data1 = load('mat/data_001.mat', 't', 'x1', 'x2');
data2 = load('mat/data_002.mat', 't', 'x1', 'x2');
figure;
hold on;
plot(data1.t, data1.x1-data1.x2);
plot(data2.t, data2.x1-data2.x2);
hold off
xlabel('Time [s]');
ylabel('Voltage [V]');
legend({'Slip-ring OFF', 'Slip-ring ON'});

3
slip-ring-test/figs/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
*.tex
*.pdf
*.svg

Binary file not shown.

After

Width:  |  Height:  |  Size: 64 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 75 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 158 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 172 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 144 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 37 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 168 KiB

BIN
slip-ring-test/figs/lpf.png Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.5 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 90 KiB

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 29 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 324 KiB

File diff suppressed because it is too large Load Diff

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 33 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 572 KiB

File diff suppressed because it is too large Load Diff

Binary file not shown.

After

Width:  |  Height:  |  Size: 35 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 166 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 217 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 36 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 44 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 149 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 158 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 44 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 45 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 190 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 60 KiB

1
slip-ring-test/img/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
*.mp4

Binary file not shown.

After

Width:  |  Height:  |  Size: 2.8 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.5 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.2 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.9 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 4.0 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 3.1 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 11 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 18 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 15 MiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.3 MiB

View File

@ -0,0 +1,3 @@
invalid color string: rgba(50, 48, 47)
(termite:25037): GLib-WARNING **: 14:10:26.520: GChildWatchSource: Exit status of a child process was requested but ECHILD was received by waitpid(). See the documentation of g_child_watch_source_new() for possible causes.

BIN
slip-ring-test/index.html Normal file

Binary file not shown.

1181
slip-ring-test/index.org Normal file

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,101 @@
% 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.
data = load('mat/data_018.mat', 'data'); data = data.data;
% Transfer function of the LPF
% We compute the transfer function from the signal without the LPF to the signal measured with the LPF.
dt = data(2, 3)-data(1, 3);
Fs = 1/dt; % [Hz]
win = hanning(ceil(10*Fs));
[Glpf, f] = tfestimate(data(:, 2), data(:, 1), win, [], [], Fs);
% We compare this transfer function with a transfer function corresponding to an ideal first order LPF with a cut-off frequency of $1000rad/s$.
% We obtain the result on figure [[fig:Glpf_bode]].
Gth = 1/(1+s/1000)
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(Glpf));
plot(f, abs(squeeze(freqresp(Gth, f, 'Hz'))));
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude');
ax2 = subplot(2, 1, 2);
hold on;
plot(f, mod(180+180/pi*phase(Glpf), 360)-180);
plot(f, 180/pi*unwrap(angle(squeeze(freqresp(Gth, f, 'Hz')))));
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]);
% Load data
% We load the data of the z axis of two geophones.
data = load('mat/data_019.mat', 'data'); data = data.data;
% Transfer function of the LPF
% We compute the transfer function from the signal without the LPF to the signal measured with the LPF.
dt = data(2, 3)-data(1, 3);
Fs = 1/dt; % [Hz]
win = hanning(ceil(10*Fs));
[Glpf, f] = tfestimate(data(:, 2), data(:, 1), win, [], [], Fs);
% We compare this transfer function with a transfer function corresponding to an ideal first order LPF with a cut-off frequency of $1060Hz$.
% We obtain the result on figure [[fig:Glpf_bode_bis]].
Gth = 1/(1+s/1060/2/pi);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(Glpf));
plot(f, abs(squeeze(freqresp(Gth, f, 'Hz'))));
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude');
ax2 = subplot(2, 1, 2);
hold on;
plot(f, mod(180+180/pi*phase(Glpf), 360)-180);
plot(f, 180/pi*unwrap(angle(squeeze(freqresp(Gth, f, 'Hz')))));
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]);

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Some files were not shown because too many files have changed in this diff Show More