diff --git a/ground-motion/data/ground_meas_id31.zip b/ground-motion/data/ground_meas_id31.zip index 8451f09..ad8f547 100644 Binary files a/ground-motion/data/ground_meas_id31.zip and b/ground-motion/data/ground_meas_id31.zip differ diff --git a/ground-motion/figs/geophone_sensibility.png b/ground-motion/figs/geophone_sensibility.png index 6748a7e..048221f 100644 Binary files a/ground-motion/figs/geophone_sensibility.png and b/ground-motion/figs/geophone_sensibility.png differ diff --git a/ground-motion/figs/ground_motion_id31_psd_displacement.png b/ground-motion/figs/ground_motion_id31_psd_displacement.png index a3129b3..879d8cd 100644 Binary files a/ground-motion/figs/ground_motion_id31_psd_displacement.png and b/ground-motion/figs/ground_motion_id31_psd_displacement.png differ diff --git a/ground-motion/figs/time_domain_displacement.png b/ground-motion/figs/time_domain_displacement.png new file mode 100644 index 0000000..8941372 Binary files /dev/null and b/ground-motion/figs/time_domain_displacement.png differ diff --git a/ground-motion/figs/time_domain_velocity.png b/ground-motion/figs/time_domain_velocity.png new file mode 100644 index 0000000..64f59b6 Binary files /dev/null and b/ground-motion/figs/time_domain_velocity.png differ diff --git a/ground-motion/index.html b/ground-motion/index.html index 17276bc..039e221 100644 Binary files a/ground-motion/index.html and b/ground-motion/index.html differ diff --git a/ground-motion/index.org b/ground-motion/index.org index a795cd0..585b73f 100644 --- a/ground-motion/index.org +++ b/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) -<> + <> #+end_src #+begin_src matlab :exports none :results silent :noweb yes -<> + <> #+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") +<> +#+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") +<> +#+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$]'); diff --git a/ground-motion/mat/psd_gm.mat b/ground-motion/mat/psd_gm.mat new file mode 100644 index 0000000..ea78259 Binary files /dev/null and b/ground-motion/mat/psd_gm.mat differ diff --git a/ground-motion/matlab/ground_meas_id31.m b/ground-motion/matlab/ground_meas_id31.m index 0a9ec39..87cd5b2 100644 --- a/ground-motion/matlab/ground_meas_id31.m +++ b/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$]');