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