Update ground motion analysis

ground-motion/index.org

@@ -74,11 +74,11 @@ On figure [[fig:ground_motion_measurements]] is an overview of multiple measurem
 
 ** 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>>
+  <<matlab-dir>>
 #+end_src
 
 #+begin_src matlab :exports none :results silent :noweb yes
-<<matlab-init>>
+  <<matlab-init>>
 #+end_src
 
 ** Load data
@@ -86,7 +86,7 @@ On figure [[fig:ground_motion_measurements]] is an overview of multiple measurem
   data = load('mat/data_028.mat', 'data'); data = data.data;
 #+end_src
 
-** Time domain plots
+** Time domain plots of the measured voltage
 #+begin_src matlab
   figure;
   hold on;
@@ -151,10 +151,16 @@ The Geophone used are L22. Their sensibility is shown on figure [[fig:geophone_s
   S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
 #+end_src
 
-#+begin_src matlab :results none :exports none
+#+begin_src matlab :results none
+  freqs = logspace(-1, 2, 1000);
+
   figure;
-  bodeFig({S}, logspace(-1, 2, 1000));
-  ylabel('Amplitude $\left[\frac{V}{m/s}\right]$')
+  hold on;
+  plot(f, abs(squeeze(freqresp(S, f, 'Hz'))));
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Magnitude $\left[\frac{V}{m/s}\right]$');
+  xlim([0.1, 100]);
 #+end_src
 
 #+NAME: fig:geophone_sensibility
@@ -176,20 +182,70 @@ We also take into account the gain of the electronics which is here set to be $6
   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}$.
+We divide the PSD measured (in $\text{V^2}/\sqrt{Hz}$) by the square of the gain of the voltage amplifier to obtain the PSD of the voltage across the geophone.
+We further divide the result by the square of the magnitude of sensibility of the Geophone to obtain the PSD of the velocity in $(m/s)^2/Hz$.
 
 #+begin_src matlab :results none
-  scaling = 1./squeeze(abs(freqresp(G0*S, f, 'Hz')));
+  psd_gv = px_dc./abs(squeeze(freqresp(G0*S, f, 'Hz'))).^2;
 #+end_src
 
+Finally, we obtain the PSD of the ground motion in $m^2/Hz$ by dividing by the square of the frequency in $rad/s$.
+#+begin_src matlab
+  psd_gm = psd_gv./(2*pi*f).^2;
+#+end_src
+
+** Time domain plots of the ground motion
+We can inverse the dynamics of the geophone to convert the measured voltage into the estimated ground motion.
+
+#+begin_src matlab
+  est_vel = lsim(inv(G0*S)*(s/2/pi)/(1+s/2/pi), data(:, 1), data(:, 3)); % Estimated velocity above 1Hz
+  est_vel = est_vel - mean(est_vel(data(:,3)>10)); % The mean value of the velocity if removed
+  est_dis = lsim(1/(1+s/2/pi), est_vel, data(:, 3)); % The velocity is integrated above 1Hz
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(data(:, 3), est_vel);
+  hold off;
+  xlabel('Time [s]'); ylabel('Velocity [m/s]');
+  xlim([10, 100]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/time_domain_velocity.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+<<plt-matlab>>
+#+end_src
+
+#+NAME: fig:time_domain_velocity
+#+CAPTION: Time domain velocity ([[./figs/time_domain_velocity.png][png]], [[./figs/time_domain_velocity.pdf][pdf]])
+[[file:figs/time_domain_velocity.png]]
+
+#+begin_src matlab :exports none
+  figure;
+  hold on;
+  plot(data(:, 3), est_dis);
+  hold off;
+  xlabel('Time [s]'); ylabel('Displacement [m]');
+  xlim([10, 100]);
+#+end_src
+
+#+HEADER: :tangle no :exports results :results none :noweb yes
+#+begin_src matlab :var filepath="figs/time_domain_displacement.pdf" :var figsize="wide-normal" :post pdf2svg(file=*this*, ext="png")
+<<plt-matlab>>
+#+end_src
+
+#+NAME: fig:time_domain_displacement
+#+CAPTION: Time domain absolute displacement ([[./figs/time_domain_displacement.png][png]], [[./figs/time_domain_displacement.pdf][pdf]])
+[[file:figs/time_domain_displacement.png]]
+
 ** Computation of the ASD of the velocity and displacement
 The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd_velocity]].
 
 #+begin_src matlab :results none
   figure;
   hold on;
-  plot(f, sqrt(px_dc).*scaling);
+  plot(f, sqrt(psd_gv));
   hold off;
   set(gca, 'xscale', 'log');
   set(gca, 'yscale', 'log');
@@ -209,11 +265,10 @@ The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd
 [[file:figs/ground_motion_id31_asd_velocity.png]]
 
 We also plot the ASD in displacement (figure [[fig:ground_motion_id31_asd_displacement]]);
-
 #+begin_src matlab :results none
   figure;
   hold on;
-  plot(f, (sqrt(px_dc).*scaling)./(2*pi*f));
+  plot(f, sqrt(psd_gm));
   hold off;
   set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
   xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
@@ -237,12 +292,12 @@ One can then compare this curve with the figure [[fig:ground_motion_measurements
 #+begin_src matlab :results none
   figure;
   hold on;
-  plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f).*1e6).^2);
+  plot(f, psd_gm.*1e12);
   hold off;
   set(gca, 'xscale', 'log');
   set(gca, 'yscale', 'log');
   xlabel('Frequency [Hz]'); ylabel('PSD of the measured displacement $\left[\frac{{ \mu m }^2}{Hz}\right]$')
-  xlim([0.1, 500]);
+  xlim([0.1, 500]); ylim([1e-13, 1e3]);
 #+end_src
 
 #+NAME: fig:ground_motion_id31_psd_displacement
@@ -256,6 +311,12 @@ One can then compare this curve with the figure [[fig:ground_motion_measurements
 #+RESULTS: fig:ground_motion_id31_psd_displacement
 [[file:figs/ground_motion_id31_psd_displacement.png]]
 
+** Save
+We save the PSD of the ground motion for further analysis.
+#+begin_src matlab
+  save('./mat/psd_gm.mat', 'f', 'psd_gm', 'psd_gv');
+#+end_src
+
 ** Comparison with other measurements of ground motion
 Now we will compare with other measurements.
 
@@ -289,8 +350,7 @@ And we compare all the measurements (figure [[fig:ground_motion_compare]]).
   hold on;
   plot(id09_f, id09_pxx, 'DisplayName', 'ID09');
   plot(cern_f, cern_pxx, 'DisplayName', 'CERN');
-  plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f)).^2, 'k', 'DisplayName', 'ID31');
-  plot(f, px_disp);
+  plot(f, psd_gm, 'k', 'DisplayName', 'ID31');
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   xlabel('Frequency [Hz]'); ylabel('PSD [$m^2/Hz$]');

ground-motion/matlab/ground_meas_id31.m

@@ -44,9 +44,15 @@ f0 = 2; % Cut-off frequency [Hz]
 
 S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
 
+freqs = logspace(-1, 2, 1000);
+
 figure;
-bodeFig({S}, logspace(-1, 2, 1000));
-ylabel('Amplitude $\left[\frac{V}{m/s}\right]$')
+hold on;
+plot(f, abs(squeeze(freqresp(S, f, 'Hz'))));
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude $\left[\frac{V}{m/s}\right]$');
+xlim([0.1, 100]);
 
 
 
@@ -64,19 +70,25 @@ G0 = 10^(G0_db/20); % [abs]
 
 
 
-% 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}$.
+% We divide the PSD measured (in $\text{V^2}/\sqrt{Hz}$) by the square of the gain of the voltage amplifier to obtain the PSD of the voltage across the geophone.
+% We further divide the result by the square of the magnitude of sensibility of the Geophone to obtain the PSD of the velocity in $(m/s)^2/Hz$.
 
 
-scaling = 1./squeeze(abs(freqresp(G0*S, f, 'Hz')));
+psd_gv = px_dc./abs(squeeze(freqresp(G0*S, f, 'Hz'))).^2;
 
-% Computation of the ASD of the velocity
+
+
+% Finally, we obtain the PSD of the ground motion in $m^2/Hz$ by dividing by the square of the frequency in $rad/s$.
+
+psd_gm = psd_gv./(2*pi*f).^2;
+
+% Computation of the ASD of the velocity and displacement
 % The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd_velocity]].
 
 
 figure;
 hold on;
-plot(f, sqrt(px_dc).*scaling);
+plot(f, sqrt(psd_gv));
 hold off;
 set(gca, 'xscale', 'log');
 set(gca, 'yscale', 'log');
@@ -92,10 +104,9 @@ xlim([0.1, 500]);
 
 % We also plot the ASD in displacement (figure [[fig:ground_motion_id31_asd_displacement]]);
 
-
 figure;
 hold on;
-plot(f, (sqrt(px_dc).*scaling)./(2*pi*f));
+plot(f, sqrt(psd_gm));
 hold off;
 set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
 xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
@@ -107,16 +118,24 @@ xlim([0.1, 500]);
 % #+CAPTION: Amplitude Spectral Density of the Displacement
 % #+RESULTS: fig:ground_motion_id31_asd_displacement
 % [[file:figs/ground_motion_id31_asd_displacement.png]]
-% And also in $\frac{{\mu u}^2}{Hz}$ (figure [[fig:ground_motion_id31_psd_displacement]]).
+
+% And we also plot the PSD of the displacement in $\frac{{\mu u}^2}{Hz}$ as it is a usual unit used (figure [[fig:ground_motion_id31_psd_displacement]]).
+% One can then compare this curve with the figure [[fig:ground_motion_measurements]].
+
 
 figure;
 hold on;
-plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f).*1e6).^2);
+plot(f, psd_gm.*1e12);
 hold off;
 set(gca, 'xscale', 'log');
 set(gca, 'yscale', 'log');
 xlabel('Frequency [Hz]'); ylabel('PSD of the measured displacement $\left[\frac{{ \mu m }^2}{Hz}\right]$')
-xlim([0.1, 500]);
+xlim([0.1, 500]); ylim([1e-13, 1e3]);
+
+% Save
+% We save the PSD of the ground motion for further analysis.
+
+save('./mat/psd_gm', 'f', 'psd_gm');
 
 % Load the measurement data
 % First we load the measurement data.
@@ -144,7 +163,7 @@ figure;
 hold on;
 plot(id09_f, id09_pxx, 'DisplayName', 'ID09');
 plot(cern_f, cern_pxx, 'DisplayName', 'CERN');
-plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f)).^2, 'k', 'DisplayName', 'ID31');
+plot(f, psd_gm, 'k', 'DisplayName', 'ID31');
 hold off;
 set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 xlabel('Frequency [Hz]'); ylabel('PSD [$m^2/Hz$]');