154 lines
4.3 KiB
Matlab
154 lines
4.3 KiB
Matlab
%% Clear Workspace and Close figures
|
|
clear; close all; clc;
|
|
|
|
%% Intialize Laplace variable
|
|
s = zpk('s');
|
|
|
|
% Load data
|
|
|
|
m_z = load('mat/data_037.mat', 'data'); m_z = m_z.data;
|
|
m_n = load('mat/data_038.mat', 'data'); m_n = m_n.data;
|
|
m_e = load('mat/data_039.mat', 'data'); m_e = m_e.data;
|
|
|
|
% Time domain plots
|
|
|
|
figure;
|
|
subplot(1, 3, 1);
|
|
hold on;
|
|
plot(m_z(:, 3), m_z(:, 2), 'DisplayName', 'Marble - Z');
|
|
plot(m_z(:, 3), m_z(:, 1), 'DisplayName', 'Floor - Z');
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Voltage [V]');
|
|
xlim([0, 100]); ylim([-2 2]);
|
|
legend('Location', 'northeast');
|
|
|
|
subplot(1, 3, 2);
|
|
hold on;
|
|
plot(m_n(:, 3), m_n(:, 2), 'DisplayName', 'Marble - N');
|
|
plot(m_n(:, 3), m_n(:, 1), 'DisplayName', 'Floor - N');
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Voltage [V]');
|
|
xlim([0, 100]); ylim([-2 2]);
|
|
legend('Location', 'northeast');
|
|
|
|
subplot(1, 3, 3);
|
|
hold on;
|
|
plot(m_e(:, 3), m_e(:, 2), 'DisplayName', 'Marble - E');
|
|
plot(m_e(:, 3), m_e(:, 1), 'DisplayName', 'Floor - E');
|
|
hold off;
|
|
xlabel('Time [s]'); ylabel('Voltage [V]');
|
|
xlim([0, 100]); ylim([-2 2]);
|
|
legend('Location', 'northeast');
|
|
|
|
% Compute the power spectral densities
|
|
% We first compute some parameters that will be used for the PSD computation.
|
|
|
|
dt = m_z(2, 3)-m_z(1, 3);
|
|
|
|
Fs = 1/dt; % [Hz]
|
|
|
|
win = hanning(ceil(10*Fs));
|
|
|
|
|
|
|
|
% Then we compute the Power Spectral Density using =pwelch= function.
|
|
|
|
[px_fz, f] = pwelch(m_z(:, 1), win, [], [], Fs);
|
|
[px_gz, ~] = pwelch(m_z(:, 2), win, [], [], Fs);
|
|
|
|
[px_fn, ~] = pwelch(m_n(:, 1), win, [], [], Fs);
|
|
[px_gn, ~] = pwelch(m_n(:, 2), win, [], [], Fs);
|
|
|
|
[px_fe, ~] = pwelch(m_e(:, 1), win, [], [], Fs);
|
|
[px_ge, ~] = pwelch(m_e(:, 2), win, [], [], Fs);
|
|
|
|
|
|
|
|
|
|
% The results are shown on figure [[fig:floor_marble_psd_z]] for the Z direction, figure [[fig:floor_marble_psd_n]] for the north direction, and figure [[fig:floor_marble_psd_e]] for the east direction.
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(px_fz), 'DisplayName', 'Floor - Z');
|
|
plot(f, sqrt(px_gz), 'DisplayName', 'Granite - Z');
|
|
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]$')
|
|
legend('Location', 'southwest');
|
|
xlim([0.1, 500]);
|
|
|
|
|
|
|
|
% #+NAME: fig:floor_marble_psd_z
|
|
% #+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - Z direction
|
|
% #+RESULTS: fig:floor_marble_psd_z
|
|
% [[file:figs/floor_marble_psd_z.png]]
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(px_fn), 'DisplayName', 'Floor - N');
|
|
plot(f, sqrt(px_gn), 'DisplayName', 'Granite - N');
|
|
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]$')
|
|
legend('Location', 'southwest');
|
|
xlim([0.1, 500]);
|
|
|
|
|
|
|
|
% #+NAME: fig:floor_marble_psd_n
|
|
% #+CAPTION: Amplitude Spectral Density of the measured voltage corresponding to the geophone located on the floor and on the marble - N direction
|
|
% #+RESULTS: fig:floor_marble_psd_n
|
|
% [[file:figs/floor_marble_psd_n.png]]
|
|
|
|
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(px_fe), 'DisplayName', 'Floor - E');
|
|
plot(f, sqrt(px_ge), 'DisplayName', 'Granite - E');
|
|
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]$')
|
|
legend('Location', 'southwest');
|
|
xlim([0.1, 500]);
|
|
|
|
% Compute the transfer function from floor motion to ground motion
|
|
% We now compute the transfer function from the floor motion to the granite motion.
|
|
|
|
% The result is shown on figure [[fig:tf_granite]].
|
|
|
|
[TZ, f] = tfestimate(m_z(:, 1), -m_z(:, 2), win, [], [], Fs);
|
|
[TN, ~] = tfestimate(m_n(:, 1), -m_n(:, 2), win, [], [], Fs);
|
|
[TE, ~] = tfestimate(m_e(:, 1), -m_e(:, 2), win, [], [], Fs);
|
|
|
|
figure;
|
|
ax1 = subplot(2, 1, 1);
|
|
hold on;
|
|
plot(f, abs(TZ), 'DisplayName', 'Z');
|
|
plot(f, abs(TN), 'DisplayName', 'N');
|
|
plot(f, abs(TE), 'DisplayName', 'E');
|
|
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(f, mod(180+180/pi*phase(TZ), 360)-180);
|
|
plot(f, mod(180+180/pi*phase(TN), 360)-180);
|
|
plot(f, mod(180+180/pi*phase(TE), 360)-180);
|
|
hold off;
|
|
set(gca, 'xscale', 'log');
|
|
ylim([-180, 180]);
|
|
yticks([-180, -90, 0, 90, 180]);
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([10, 100]);
|