%% Clear Workspace and Close figures
clear; close all; clc;

%% Intialize Laplace variable
s = zpk('s');

% Load data
% We load the data of the z axis of two geophones.


meas_sr = load('mat/data_018.mat', 'data'); meas_sr = meas_sr.data;
meas_di = load('mat/data_019.mat', 'data'); meas_di = meas_di.data;

% Analysis - Time Domain
% First, we compare the time domain signals for the two experiments (figure [[fig:slipring_time]]).


figure;
hold on;
plot(meas_di(:, 3), meas_di(:, 2), 'DisplayName', 'Geophone - Direct');
plot(meas_sr(:, 3), meas_sr(:, 2), 'DisplayName', 'Geophone - Slip-Ring');
hold off;
xlabel('Time [s]'); ylabel('Voltage [V]');
xlim([0, 50]);
legend('location', 'northeast');

% Analysis - Frequency Domain
% We then compute the Power Spectral Density of the two signals and we compare them (figure [[fig:slipring_asd]]).


dt = meas_di(2, 3) - meas_di(1, 3);
Fs = 1/dt;

win = hanning(ceil(5*Fs));

[px_di, f] = pwelch(meas_di(:, 2), win, [], [], Fs);
[px_sr, ~] = pwelch(meas_sr(:, 2), win, [], [], Fs);

figure;
hold on;
plot(f, sqrt(px_sr), 'DisplayName', 'Slip-Ring');
plot(f, sqrt(px_di), 'DisplayName', 'Wire');
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density $\left[\frac{V}{\sqrt{Hz}}\right]$')
xlim([1, 500]);
legend('Location', 'southwest');