% Matlab Init                                              :noexport:ignore:

clear; close all; clc;

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

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

% Load data
% We first load the data for the three axis.

z     = load('mat/data_001.mat', 't', 'x1', 'x2');
east  = load('mat/data_002.mat', 't', 'x1', 'x2');
north = load('mat/data_003.mat', 't', 'x1', 'x2');

% Compare PSD
% The PSD for each axis of the two geophones are computed.

[pz1, fz] = pwelch(z.x1, hanning(ceil(length(z.x1)/100)), [], [], 1/(z.t(2)-z.t(1)));
[pz2, ~]  = pwelch(z.x2, hanning(ceil(length(z.x2)/100)), [], [], 1/(z.t(2)-z.t(1)));

[pe1, fe] = pwelch(east.x1, hanning(ceil(length(east.x1)/100)), [], [], 1/(east.t(2)-east.t(1)));
[pe2, ~]  = pwelch(east.x2, hanning(ceil(length(east.x2)/100)), [], [], 1/(east.t(2)-east.t(1)));

[pn1, fn] = pwelch(north.x1, hanning(ceil(length(north.x1)/100)), [], [], 1/(north.t(2)-north.t(1)));
[pn2, ~]  = pwelch(north.x2, hanning(ceil(length(north.x2)/100)), [], [], 1/(north.t(2)-north.t(1)));



% We compare them. The result is shown on figure [[fig:compare_axis_psd]].

figure;
hold on;
plot(fz, sqrt(pz1), '-',   'Color', [0 0.4470 0.7410], 'DisplayName', 'z');
plot(fz, sqrt(pz2), '--',  'Color', [0 0.4470 0.7410], 'HandleVisibility', 'off');
plot(fe, sqrt(pe1), '-',   'Color', [0.8500 0.3250 0.0980], 'DisplayName', 'east');
plot(fe, sqrt(pe2), '--',  'Color', [0.8500 0.3250 0.0980], 'HandleVisibility', 'off');
plot(fn, sqrt(pn1), '-',   'Color', [0.9290 0.6940 0.1250], 'DisplayName', 'north');
plot(fn, sqrt(pn2), '--',  'Color', [0.9290 0.6940 0.1250], 'HandleVisibility', 'off');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]');
legend('Location', 'northeast');
xlim([0, 500]);

% Compare TF
% The transfer functions from one geophone to the other are also computed for each axis.
% The result is shown on figure [[fig:compare_tf_axis]].


[Tz, fz] = tfestimate(z.x1, z.x2, hanning(ceil(length(z.x1)/100)), [], [], 1/(z.t(2)-z.t(1)));
[Te, fe] = tfestimate(east.x1, east.x2, hanning(ceil(length(east.x1)/100)), [], [], 1/(east.t(2)-east.t(1)));
[Tn, fn] = tfestimate(north.x1, north.x2, hanning(ceil(length(north.x1)/100)), [], [], 1/(north.t(2)-north.t(1)));

figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(fz, abs(Tz), 'DisplayName', 'z');
plot(fe, abs(Te), 'DisplayName', 'east');
plot(fn, abs(Tn), 'DisplayName', 'north');
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(fz, mod(180+180/pi*phase(Tz), 360)-180);
plot(fe, mod(180+180/pi*phase(Te), 360)-180);
plot(fn, mod(180+180/pi*phase(Tn), 360)-180);
hold off;
set(gca, 'xscale', 'log');
ylim([-180, 180]);
yticks([-180, -90, 0, 90, 180]);
xlabel('Frequency [Hz]'); ylabel('Phase');

linkaxes([ax1,ax2],'x');
xlim([1, 500]);