%% Clear Workspace and Close figures
clear; close all; clc;

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

% Signal's PSD
% We load the PSD of the signal we wish to replicate.

load('./mat/dist_psd.mat', 'dist_f');



% We remove the first value with very high PSD.

dist_f.f = dist_f.f(3:end);
dist_f.psd_gm = dist_f.psd_gm(3:end);



% The PSD of the signal is shown on figure ref:fig:psd_ground_motion.

figure;
hold on;
plot(dist_f.f, dist_f.psd_gm)
hold off;
xlabel('Frequency [Hz]');
ylabel('Power Spectral Density');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([0.1, 500]);

% Transfer Function that approximate the ASD
% Using =sisotool= or any other tool, we create a transfer function $G$ such that its magnitude is close to the Amplitude Spectral Density $\Gamma_x = \sqrt{S_x}$:
% \[ |G(j\omega)| \approx \Gamma_x(\omega) \]


G_gm = 0.002*(s^2 + 3.169*s + 27.74)/(s*(s+32.73)*(s+8.829)*(s+7.983)^2);



% We compare the ASD $\Gamma_x(\omega)$ and the magnitude of the generated transfer function $|G(j\omega)|$ in figure [[]].

figure;
hold on;
plot(dist_f.f, sqrt(dist_f.psd_gm))
plot(dist_f.f, abs(squeeze(freqresp(G_gm, dist_f.f, 'Hz'))))
hold off;
xlabel('Frequency [Hz]');
ylabel('Amplitude Spectral Density');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([0.1, 500]);

% Generated Time domain signal
% We know that a signal $u$ going through a LTI system $G$ (figure [[fig:velocity_to_voltage_psd_bis]]) will have its ASD modified according to the following equation:
% \begin{equation}
%   \Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{uu}(\omega)
% \end{equation}

% #+NAME: fig:velocity_to_voltage_psd_bis
% #+CAPTION: Schematic of the instrumentation used for the measurement
% [[file:figs/velocity_to_voltage_psd.png]]

% Thus, if we create a random signal with an ASD equal to one at all frequency and we pass this signal through the previously defined LTI transfer function =G_gm=, we should obtain a signal with an ASD that approximate the original ASD.

% To obtain a random signal with an ASD equal to one, we use the following code.

Fs = 2*dist_f.f(end); % Sampling Frequency [Hz]
Ts = 1/Fs; % Sampling Time [s]

t = 0:Ts:500; % Time Vector [s]
u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one



% We then use =lsim= to compute $y$ as shown in figure [[fig:velocity_to_voltage_psd_bis]].

y = lsim(G_gm, u, t);



% The obtained time domain signal is shown in figure [[fig:time_domain_u]].

figure;
plot(t, y);
xlabel('Time [s]');
ylabel('Amplitude');

% Comparison of the Power Spectral Densities
% We now compute the Power Spectral Density of the computed time domain signal $y$.

nx = length(y);
na = 16;
win = hanning(floor(nx/na));

[pxx, f] = pwelch(y, win, 0, [], Fs);



% Finally, we compare the PSD of the original signal and the obtained signal on figure ref:fig:psd_comparison.

figure;
hold on;
plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
plot(f, pxx, 'DisplayName', 'Computed')
hold off;
xlabel('Frequency [Hz]');
ylabel('Power Spectral Density');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
legend('location', 'northeast');
xlim([0.1, 500]);

% Simulink
% One advantage of this technique is that it can be easily integrated into simulink.

% The corresponding schematic is shown in figure [[fig:simulink_psd_generate]] where the block =Band-Limited White Noise= is used to generate a random signal with a PSD equal to one (parameter =Noise Power= is set to 1).

% Then, the signal generated pass through the transfer function representing the wanted ASD.

% #+name: fig:simulink_psd_generate
% #+caption: Simulink Schematic
% [[file:figs/simulink_psd_generate.png]]

% We simulate the system shown in figure [[fig:simulink_psd_generate]].

out = sim('matlab/generate_signal_psd.slx');



% And we compute the PSD of the generated signal.

nx = length(out.u_gm.Data);
na = 8;
win = hanning(floor(nx/na));

[pxx, f] = pwelch(out.u_gm.Data, win, 0, [], 1e3);



% Finally, we compare the PSD of the generated signal with the original PSD in figure [[fig:compare_psd_original_simulink]].

figure;
hold on;
plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
plot(f, pxx, 'DisplayName', 'Computed from Simulink')
hold off;
xlabel('Frequency [Hz]');
ylabel('Power Spectral Density');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
legend('location', 'northeast');
xlim([0.1, 500]);