% Matlab Init                                             :noexport:ignore:

clear; close all; clc;

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

%% Initialize ans with org-babel
ans = 0;

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

sr_off = load('mat/data_012.mat', 'data'); sr_off = sr_off.data;
sr_on  = load('mat/data_013.mat', 'data'); sr_on  = sr_on.data;

% Time Domain
% We compare the signal when the Slip-Ring is OFF (figure [[fig:sr_geophone_time_off]]) and when it is ON (figure [[fig:sr_geophone_time_on]]).


figure;
hold on;
plot(sr_off(:, 3), sr_off(:, 1), 'DisplayName', 'Direct');
plot(sr_off(:, 3), sr_off(:, 2), 'DisplayName', 'Slip-Ring');
hold off;
legend('Location', 'northeast');
xlabel('Time [s]');
ylabel('Voltage [V]');



% #+NAME: fig:sr_geophone_time_off
% #+CAPTION: Comparison of the time domain signals when the slip-ring is OFF
% #+RESULTS: fig:sr_geophone_time_off
% [[file:figs/sr_geophone_time_off.png]]


figure;
hold on;
plot(sr_on(:, 3),  sr_on(:, 1),  'DisplayName', 'Direct');
plot(sr_on(:, 3),  sr_on(:, 2),  'DisplayName', 'Slip-Ring');
hold off;
legend('Location', 'northeast');
xlabel('Time [s]');
ylabel('Voltage [V]');

% Frequency Domain
% We first compute some parameters that will be used for the PSD computation.

dt = sr_off(2, 3)-sr_off(1, 3);

Fs = 1/dt; % [Hz]

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



% Then we compute the Power Spectral Density using =pwelch= function.

% Direct measure
[pxdoff, ~] = pwelch(sr_off(:, 1), win, [], [], Fs);
[pxdon,  ~] = pwelch(sr_on(:, 1),  win, [], [], Fs);

% Slip-Ring measure
[pxsroff, f] = pwelch(sr_off(:, 2), win, [], [], Fs);
[pxsron,  ~] = pwelch(sr_on(:, 2),  win, [], [], Fs);



% Finally, we compare the Amplitude Spectral Density of the signals (figure [[fig:sr_geophone_asd]]);


figure;
hold on;
plot(f, sqrt(pxdoff), 'DisplayName', 'Direct - OFF');
plot(f, sqrt(pxsroff), 'DisplayName', 'Slip-Ring - OFF');
plot(f, sqrt(pxdon),  'DisplayName', 'Direct - ON');
plot(f, sqrt(pxsron),  'DisplayName', 'Slip-Ring - ON');
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', 'northeast');
xlim([0.1, 500]);



% #+NAME: fig:sr_geophone_asd
% #+CAPTION: Comparison of the Amplitude Spectral Sensity
% #+RESULTS: fig:sr_geophone_asd
% [[file:figs/sr_geophone_asd.png]]


xlim([100, 500]);

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

sr_lpf_off = load('mat/data_016.mat', 'data'); sr_lpf_off = sr_lpf_off.data;
sr_lpf_on  = load('mat/data_017.mat', 'data'); sr_lpf_on  = sr_lpf_on.data;

% Time Domain
% We compare the signal when the Slip-Ring is OFF (figure [[fig:sr_lpf_geophone_time_off]]) and when it is ON (figure [[fig:sr_lpf_geophone_time_on]]).


figure;
hold on;
plot(sr_lpf_off(:, 3), sr_lpf_off(:, 1), 'DisplayName', 'Direct');
plot(sr_lpf_off(:, 3), sr_lpf_off(:, 2), 'DisplayName', 'Slip-Ring');
hold off;
legend('Location', 'northeast');
xlabel('Time [s]');
ylabel('Voltage [V]');



% #+NAME: fig:sr_lpf_geophone_time_off
% #+CAPTION: Comparison of the time domain signals when the slip-ring is OFF
% #+RESULTS: fig:sr_lpf_geophone_time_off
% [[file:figs/sr_lpf_geophone_time_off.png]]


figure;
hold on;
plot(sr_lpf_on(:, 3),  sr_lpf_on(:, 1),  'DisplayName', 'Direct');
plot(sr_lpf_on(:, 3),  sr_lpf_on(:, 2),  'DisplayName', 'Slip-Ring');
hold off;
legend('Location', 'northeast');
xlabel('Time [s]');
ylabel('Voltage [V]');

% Frequency Domain
% We first compute some parameters that will be used for the PSD computation.

dt = sr_lpf_off(2, 3)-sr_lpf_off(1, 3);

Fs = 1/dt; % [Hz]

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



% Then we compute the Power Spectral Density using =pwelch= function.

% Direct measure
[pxd_lpf_off, ~] = pwelch(sr_lpf_off(:, 1), win, [], [], Fs);
[pxd_lpf_on,  ~] = pwelch(sr_lpf_on(:, 1),  win, [], [], Fs);

% Slip-Ring measure
[pxsr_lpf_off, f] = pwelch(sr_lpf_off(:, 2), win, [], [], Fs);
[pxsr_lpf_on,  ~] = pwelch(sr_lpf_on(:, 2),  win, [], [], Fs);



% Finally, we compare the Amplitude Spectral Density of the signals (figure [[fig:sr_lpf_geophone_asd]]);


figure;
hold on;
plot(f, sqrt(pxd_lpf_off), 'DisplayName', 'Direct - OFF');
plot(f, sqrt(pxsr_lpf_off), 'DisplayName', 'Slip-Ring - OFF');
plot(f, sqrt(pxd_lpf_on),  'DisplayName', 'Direct - ON');
plot(f, sqrt(pxsr_lpf_on),  'DisplayName', 'Slip-Ring - ON');
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', 'northeast');
xlim([0.1, 500]);



% #+NAME: fig:sr_lpf_geophone_asd
% #+CAPTION: Comparison of the Amplitude Spectral Sensity
% #+RESULTS: fig:sr_lpf_geophone_asd
% [[file:figs/sr_lpf_geophone_asd.png]]


xlim([100, 500]);

% Comparison of with and without LPF

figure;
hold on;
plot(f, sqrt(pxdon),  'DisplayName', 'Direct - ON');
plot(f, sqrt(pxsron),  'DisplayName', 'Slip-Ring - ON');
plot(f, sqrt(pxd_lpf_on),  'DisplayName', 'Direct - ON - LPF');
plot(f, sqrt(pxsr_lpf_on),  'DisplayName', 'Slip-Ring - ON - LPF');
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', 'northeast');
xlim([0.1, 500]);