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

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

addpath('./mat/');

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

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]);

% Algorithm
% We define some parameters that will be used in the algorithm.

Fs = 2*dist_f.f(end);    % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
N  = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD
T0 = N/Fs;               % Signal Duration [s]
df = 1/T0;               % Frequency resolution of the DFT [Hz]
                         % Also equal to (dist_f.f(2)-dist_f.f(1))



% We then specify the wanted PSD.

phi = dist_f.psd_gm;



% We create amplitudes corresponding to wanted PSD.

C = zeros(N/2,1);
for i = 1:N/2
  C(i) = sqrt(phi(i)*df);
end



% Finally, we add some random phase to =C=.

theta = 2*pi*rand(N/2,1); % Generate random phase [rad]

Cx = [0 ; C.*complex(cos(theta),sin(theta))];
Cx = [Cx; flipud(conj(Cx(2:end)))];;

% Obtained Time Domain Signal
% The time domain data is generated by an inverse FFT.

% The =ifft= Matlab does not take into account the sampling frequency, thus we need to normalize the signal.

u = N/sqrt(2)*ifft(Cx); % Normalisation of the IFFT
t = linspace(0, T0, N+1); % Time Vector [s]

figure;
plot(t, u)
xlabel('Time [s]');
ylabel('Amplitude');
xlim([t(1), t(end)]);

% PSD Comparison
% We duplicate the time domain signal to have a longer signal and thus a more precise PSD result.

u_rep = repmat(u, 10, 1);



% We compute the PSD of the obtained signal with the following commands.

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

[pxx, f] = pwelch(u_rep, 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]);