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

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

addpath('./src/');

%% Simulation configuration
Fs = 1e4; % Sampling Frequency [Hz]
Ts = 1/Fs; % Sampling Time [s]

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

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

%% Security
% u_min = -1;
% u_max = 7.5;

u_min = 1.25;
u_max = 5.25;

d_min = -50e-6;
d_max =  50e-6;

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

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

% Notch filter at the resonance of the APA
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',   400, ...
    'V_mean',  0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'sweep_type',   'log', ...
    'V_exc', G_sweep*1/(1 + s/2/pi/500));

V_sweep_lf = generateSweepExc('Ts',      Ts, ...
    'f_start', 20, ...
    'f_end',   40, ...
    'V_mean',  0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'sweep_type',   'lin', ...
    'V_exc', 0.2);

%% High Frequency Shaped Noise
[b,a] = cheby1(10, 2, 2*pi*[240 2e3], 'bandpass', 's');
wL = 0.005*tf(b, a);

V_noise_hf = generateShapedNoise('Ts', Ts, ...
    'V_mean', 0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'smooth_ends', true, ...
    'V_exc', wL);

%% Low Frequency Shaped Noise
[b,a] = cheby1(10, 2, 2*pi*[20 40], 'bandpass', 's');
wL = 0.01*tf(b, a);

V_noise_lf = generateShapedNoise('Ts', 1/Fs, ...
    'V_mean', 0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'smooth_ends', true, ...
    'V_exc', wL);

%% Band Pass Shaped Noise
[b,a] = cheby1(10, 2, 2*pi*[9 1.1e3], 'bandpass', 's');
wL = 0.005*tf(b, a);

V_noise_bf = generateShapedNoise('Ts', 1/Fs, ...
    'V_mean', 0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'smooth_ends', true, ...
    'V_exc', wL);

%% Sinus excitation with increasing amplitude
V_sin = generateSinIncreasingAmpl('Ts', 1/Fs, ...
    'V_mean', 0, ...
    'sin_ampls', [0.1, 0.2, 0.4, 1, 2, 4], ...
    'sin_period', 1, ...
    'sin_num', 5, ...
    't_start', Trec_start, ...
    'smooth_ends', true);

%% Zero Voltage
V_zero = generateShapedNoise('Ts', 1/Fs, ...
    'V_mean', 0, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'smooth_ends', true, ...
    'V_exc', tf(0));

%% Offset Voltage for all APA
V_zero_off = generateShapedNoise('Ts', 1/Fs, ...
    'V_mean', 3.25, ...
    't_start', Trec_start, ...
    'exc_duration', Trec_dur, ...
    'smooth_ends', true, ...
    'V_exc', tf(0));

%% Select the excitation signal and offset voltage
V_exc = timeseries(V_noise_bf(2,:), V_noise_bf(1,:));
V_off = timeseries(V_zero_off(2,:), V_zero_off(1,:));

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

ax1 = nexttile;
hold on;
plot(V_exc.Time, squeeze(V_exc.Data+V_off.Data));
plot([V_exc.Time(1), V_exc.Time(end)], [u_min, u_min], 'k--');
plot([V_exc.Time(1), V_exc.Time(end)], [u_max, u_max], 'k--');
hold off;
xlabel('Time [s]'); ylabel('Amplitude [V]');

ax2 = nexttile;
win = hanning(floor(length(squeeze(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 data that will be loaded in the Simulink file
save('sim_data/data_sim.mat', ...
    'Fs', 'Ts', 'Tsim', 'Trec_start', 'Trec_dur', ...
    'V_exc', 'V_off', ...
    'u_min', 'u_max', ...
    'd_min', 'd_max');