%%
clear; close all; clc;

%%
load('sim_data/data_sim.mat', 'Trec_start', 'Trec_dur', 'Ts');

%% White Noise => Identification of the sensibility transfer function
s = zpk('s');

Rx_noise = generateShapedNoise('Ts', Ts, ...
    'V_mean', 0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'V_exc', 2e-6/(2*pi*0.2 + s));

R_dist = timeseries(Rx_noise(2, :), Rx_noise(1, :));

figure;
tiledlayout(1, 2, 'TileSpacing', 'Normal', 'Padding', 'None');

ax1 = nexttile;
hold on;
plot(R_dist.Time, squeeze(R_dist.Data(1, 1, :)));
xlabel('Time [s]'); ylabel('Amplitude [m,rad]');

ax2 = nexttile;
hold on;
win = hanning(floor(length(squeeze(R_dist.Data(1, 1, :)))/8));
[pxx, f] = pwelch(squeeze(R_dist.Data(1, 1, :)), win, 0, [], 1/Ts);
plot(f, sqrt(pxx))
hold off;
xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density [$m/\sqrt{Hz}$]');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([1, 1e3]); ylim([1e-10, 1e0]);

%%
save('sim_data/ref_dist.mat', 'R_dist');