s = tf('s');
addpath('src')

%%
Fs = 10e3; % Sampling Frequency [Hz]
Ts = 1/Fs; % Sampling Time [s]

Trec_start = 5;  % Start time for Recording [s]
Trec_dur   = 100; % Recording Duration [s]

Tsim = 2*Trec_start + Trec_dur; % Simulation Time [s]

%% Sweep Sine
gc = 0.1;
xi = 0.5;
wn = 2*pi*94.3;

G_sweep = 0.2*(s^2 + 2*gc*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);

V_sweep = generateSweepExc('Ts',      Ts, ...
                           'f_start', 10, ...
                           'f_end',   2e3, ...
                           'V_mean',  3.25, ...
                           't_start', Trec_start, ...
                           'exc_duration', Trec_dur, ...
                           'sweep_type',   'log', ...
                           'V_exc', G_sweep*1/(1 + s/2/pi/500));

%% Shaped Noise
V_noise = generateShapedNoise('Ts', 1/Fs, ...
                              'V_mean', 3.25, ...
                              't_start', Trec_start, ...
                              'exc_duration', Trec_dur, ...
                              'smooth_ends', true, ...
                              'V_exc', 0.00/(1 + s/2/pi/50));

%% Select the excitation signal
V_exc = timeseries(V_noise(2,:), V_noise(1,:));

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

ax1 = nexttile;
plot(V_exc.Time, squeeze(V_exc.Data));
xlabel('Time [s]'); ylabel('Amplitude [V]');

ax2 = nexttile;
win = hanning(floor(length(V_exc.Data)/8));
[pxx, f] = pwelch(squeeze(V_exc.Data), win, 0, [], Fs);
plot(f, pxx)
xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [$V^2/Hz$]');
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
xlim([1, Fs/2]); ylim([1e-10, 1e0]);

%% Save
save('./frf_data.mat', 'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', 'V_exc');