% Matlab Init                                              :noexport:ignore:

clear; close all; clc;

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

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

% Load data

amp_off = load('mat/data_003.mat', 'data'); amp_off = amp_off.data(:, [1,3]);
amp_20d = load('mat/data_004.mat', 'data'); amp_20d = amp_20d.data(:, [1,3]);
amp_40d = load('mat/data_005.mat', 'data'); amp_40d = amp_40d.data(:, [1,3]);
amp_60d = load('mat/data_006.mat', 'data'); amp_60d = amp_60d.data(:, [1,3]);

% Time Domain
% The time domain signals are shown on figure [[fig:ampli_noise_time]].


figure;
hold on;
plot(amp_off(:, 2), amp_off(:, 1), 'DisplayName', 'OFF');
plot(amp_20d(:, 2), amp_20d(:, 1), 'DisplayName', '20dB');
plot(amp_40d(:, 2), amp_40d(:, 1), 'DisplayName', '40dB');
plot(amp_60d(:, 2), amp_60d(:, 1), 'DisplayName', '60dB');
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 = amp_off(2, 2)-amp_off(1, 2);

Fs = 1/dt; % [Hz]

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



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

[pxoff, f] = pwelch(amp_off(:,1), win, [], [], Fs);
[px20d, ~] = pwelch(amp_20d(:,1), win, [], [], Fs);
[px40d, ~] = pwelch(amp_40d(:,1), win, [], [], Fs);
[px60d, ~] = pwelch(amp_60d(:,1), win, [], [], Fs);



% We compute the theoretical ADC noise.

q = 20/2^16; % quantization
Sq = q^2/12/1000; % PSD of the ADC noise



% Finally, the ASD is shown on figure [[fig:ampli_noise_psd]].

figure;
hold on;
plot(f, sqrt(pxoff), 'DisplayName', 'OFF');
plot(f, sqrt(px20d), 'DisplayName', '20dB');
plot(f, sqrt(px40d), 'DisplayName', '40dB');
plot(f, sqrt(px60d), 'DisplayName', '60dB');
plot([0.1, 500], [sqrt(Sq), sqrt(Sq)], 'k--');
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]);