% 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.


load('mat/data_001.mat', 't', 'x1', 'x2');
dt = t(2) - t(1);

% Time Domain Data

figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([t(1), t(end)]);



% #+NAME: fig:data_time_domain
% #+CAPTION: Time domain Data
% #+RESULTS: fig:data_time_domain
% [[file:figs/data_time_domain.png]]



figure;
hold on;
plot(t, x1);
plot(t, x2);
hold off;
xlabel('Time [s]');
ylabel('Voltage [V]');
xlim([0 1]);

% Computation of the ASD of the measured voltage
% We first define the parameters for the frequency domain analysis.

win = hanning(ceil(length(x1)/100));
Fs = 1/dt;

[pxx1, f] = pwelch(x1, win, [], [], Fs);
[pxx2, ~] = pwelch(x2, win, [], [], Fs);

% Scaling to take into account the sensibility of the geophone and the voltage amplifier
% The Geophone used are L22.
% Their sensibility are shown on figure [[fig:geophone_sensibility]].


S0 = 88; % Sensitivity [V/(m/s)]
f0 = 2; % Cut-off frequnecy [Hz]
S = (s/2/pi/f0)/(1+s/2/pi/f0);

figure;
bodeFig({S});
ylabel('Amplitude [V/(m/s)]')



% #+NAME: fig:geophone_sensibility
% #+CAPTION: Sensibility of the Geophone
% #+RESULTS: fig:geophone_sensibility
% [[file:figs/geophone_sensibility.png]]


% We also take into account the gain of the electronics which is here set to be $60dB$.
% The amplifiers also include a low pass filter with a cut-off frequency set at 1kHz.


G0 = 60; % [dB]

G = 10^(G0/20)/(1+s/2/pi/1000);



% We divide the ASD measured (in $\text{V}/\sqrt{\text{Hz}}$) by the transfer function of the voltage amplifier to obtain the ASD of the voltage across the geophone.
% We further divide the result by the sensibility of the Geophone to obtain the ASD of the velocity in $m/s/\sqrt{Hz}$.


scaling = 1./squeeze(abs(freqresp(G*S, f, 'Hz')));

% Computation of the ASD of the velocity
% The ASD of the measured velocity is shown on figure [[fig:psd_velocity]].


figure;
hold on;
plot(f, sqrt(pxx1).*scaling);
plot(f, sqrt(pxx2).*scaling);
hold off;
set(gca, 'xscale', 'log');
set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
xlim([2, 500]);



% #+NAME: fig:psd_velocity
% #+CAPTION: Spectral density of the velocity
% #+RESULTS: fig:psd_velocity
% [[file:figs/psd_velocity.png]]

% We also plot the ASD in displacement (figure [[fig:asd_displacement]]);


figure;
hold on;
plot(f, (pxx1.*scaling./f).^2);
plot(f, (pxx2.*scaling./f).^2);
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [m/s/sqrt(Hz)]')
xlim([2, 500]);

% Transfer function between the two geophones
% We here compute the transfer function from one geophone to the other.
% The result is shown on figure [[fig:tf_geophones]].

% We also compute the coherence between the two signals (figure [[fig:coh_geophones]]).


[T12, ~] = tfestimate(x1, x2, win, [], [], Fs);

figure;
ax1 = subplot(2, 1, 1);
plot(f, abs(T12));
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
set(gca, 'XTickLabel',[]);
ylabel('Magnitude');

ax2 = subplot(2, 1, 2);
plot(f, mod(180+180/pi*phase(T12), 360)-180);
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]);



% #+NAME: fig:tf_geophones
% #+CAPTION: Estimated transfer function between the two geophones
% #+RESULTS: fig:tf_geophones
% [[file:figs/tf_geophones.png]]


[coh12, ~] = mscohere(x1, x2, win, [], [], Fs);

figure;
plot(f, coh12);
set(gca, 'xscale', 'log');
xlabel('Frequency [Hz]'); ylabel('Coherence');
ylim([0,1]); xlim([1, 500]);

% Estimation of the sensor noise
% The technique to estimate the sensor noise is taken from cite:barzilai98_techn_measur_noise_sensor_presen.

% The coherence between signals $X$ and $Y$ is defined as follow
% \[ \gamma^2_{XY}(\omega) = \frac{|G_{XY}(\omega)|^2}{|G_{X}(\omega)| |G_{Y}(\omega)|} \]
% where $|G_X(\omega)|$ is the output Power Spectral Density (PSD) of signal $X$ and $|G_{XY}(\omega)|$ is the Cross Spectral Density (CSD) of signal $X$ and $Y$.

% The PSD and CSD are defined as follow:
% \begin{align}
%   |G_X(\omega)|    &= \frac{2}{n_d T} \sum^{n_d}_{n=1} \left| X_k(\omega, T) \right|^2 \\
%   |G_{XY}(\omega)| &= \frac{2}{n_d T} \sum^{n_d}_{n=1} [ X_k^*(\omega, T) ] [ Y_k(\omega, T) ]
% \end{align}
% where:
% - $n_d$ is the number for records averaged
% - $T$ is the length of each record
% - $X_k(\omega, T)$ is the finite Fourier transform of the kth record
% - $X_k^*(\omega, T)$ is its complex conjugate

% The =mscohere= function is compared with this formula on Appendix (section [[sec:coherence]]), it is shown that it is identical.

% Figure [[fig:huddle_test]] illustrate a block diagram model of the system used to determine the sensor noise of the geophone.

% Two geophones are mounted side by side to ensure that they are exposed by the same motion input $U$.

% Each sensor has noise $N$ and $M$.

% #+NAME: fig:huddle_test
% #+CAPTION: Huddle test block diagram
% [[file:figs/huddle-test.png]]

% We here assume that each sensor has the same magnitude of instrumental noise ($N = M$).
% We also assume that $H_1 = H_2 = 1$.

% We then obtain:
% #+NAME: eq:coh_bis
% \begin{equation}
%   \gamma_{XY}^2(\omega) = \frac{1}{1 + 2 \left( \frac{|G_N(\omega)|}{|G_U(\omega)|} \right) + \left( \frac{|G_N(\omega)|}{|G_U(\omega)|} \right)^2}
% \end{equation}

% Since the input signal $U$ and the instrumental noise $N$ are incoherent:
% #+NAME: eq:incoherent_noise
% \begin{equation}
%   |G_X(\omega)| = |G_N(\omega)| + |G_U(\omega)|
% \end{equation}

% From equations [[eq:coh_bis]] and [[eq:incoherent_noise]], we finally obtain
% #+begin_important
% #+NAME: eq:noise_psd
% \begin{equation}
%   |G_N(\omega)| = |G_X(\omega)| \left( 1 - \sqrt{\gamma_{XY}^2(\omega)} \right)
% \end{equation}
% #+end_important

% The instrumental noise is computed below. The result in V^2/Hz is shown on figure [[fig:intrumental_noise_V]].

pxxN = pxx1.*(1 - coh12);

figure;
hold on;
plot(f, pxx1, '-');
plot(f, pxx2, '-');
plot(f, pxxN, 'k--');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [$V^2/Hz$]');
xlim([1, 500]);



% #+NAME: fig:intrumental_noise_V
% #+CAPTION: Instrumental Noise and Measurement in $V^2/Hz$
% #+RESULTS: fig:intrumental_noise_V
% [[file:figs/intrumental_noise_V.png]]

% This is then further converted into velocity and compared with the ground velocity measurement. (figure [[fig:intrumental_noise_velocity]])

figure;
hold on;
plot(f, sqrt(pxx1).*scaling, '-');
plot(f, sqrt(pxx2).*scaling, '-');
plot(f, sqrt(pxxN).*scaling, 'k--');
hold off;
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlabel('Frequency [Hz]'); ylabel('PSD [$m/s/\sqrt{Hz}$]');
xlim([1, 500]);