close all hidden
clear
clc

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 
cmap = lines;

%**************************************************************************
% Basic parameters
%**************************************************************************
f_s = 200;        % Sampling freq for "continuous" signal in Hz
T_0 = 2;          % Fundamental period for signal description in s
T_max = 100;      % Signal duration in seconds

%**************************************************************************
% Calculate support variables
%**************************************************************************
N = T_max * f_s + 1;                     % Duration in samples
t = linspace(-T_max/2, T_max/2, N);      % Time scale in s
w = 2 * pi * linspace(-f_s/2, f_s/2, N); % Frequency scale in s^-1

%**************************************************************************
% Signal generation
%**************************************************************************

signalSel = 5;

switch(signalSel)
    case 1
    % Rect signal -----------------------------------------------
    % Set time signal using ones and zeros
    v = ones(N, 1);                 % Fill signal with ones
    v(t < -T_0) = 0 * v(t < -T_0);  % Set signal before rect pulse to 1
    v(t > T_0) = 0 * v(t > T_0);    % Set signal after rect pulse to 0

    % Set discretized Fourier spectrum to sin(x)/x function
    % NaN value at zero is okay here
    V = abs(2 * sin(w * T_0) ./ w); 
    
   % Use sinc (si) function as better alternative 
   % Only with dsp toolbox
   %    V = abs(2 * T_0 * sinc(w / pi * T_0));
    
    % Calculate approx abs Fourier spectrum based on DFT
    V_dft = abs(fftshift(fft(v) / f_s));
    
    case 2
    % Sinc signal ----------------------------------------------
    % Define auxiliary parameters
    f_0 = 1;
    w_0 = 2 * pi * f_0;
    
    % Set time signal to sin(x)/x function
    % NaN value at zero is not okay here as
    % it would make the FFT result NaNs only
    v = (sin(w_0 * t)  + eps) ./ (w_0 * t + eps);
    
    % Use sinc (si) function as better alternative 
    % Only with dsp toolbox
    %  v = sinc(w_0 / pi * t);
    
    % Set discretized Fourier spectrum to rect pulse
    V = pi / w_0 * ones(N, 1);
    V(w < -w_0) = 0 * V(w < -w_0);
    V(w > w_0) = 0 * V(w > w_0);
   
    % Calculate approx abs Fourier spectrum based on DFT
    V_dft = abs(fftshift(fft(v) / f_s));
    
    case 3
    % Exp signal ----------------------------------------------
    % Define auxiliary parameters
    alpha = 0.1;
    
    % Set time signal to exp function for t > 0
    v = zeros(N, 1);
    v(t >= 0) = exp(-alpha * t(t >= 0));
    
    % Set discretized Fourier spectrum to 1/(a + jw)
    V = abs(1 ./ (alpha + 1i * w)); 
   
    % Calculate approx abs Fourier spectrum based on DFT
    V_dft = abs(fftshift(fft(v) / f_s));
    
    case 4
    % Cos signal ----------------------------------------------
    % Define auxiliary parameters
    f_0 = 1;
    w_0 = 2 * pi * f_0;
    
    % Set time signal to exp function for t > 0
    v = cos(w_0 * t);
    
    % Set discretized Fourier spectrum to zero 
    V = zeros(N, 1);
    % Find freq index closest to w_0
    [~, w_0_c] = min(abs(w - w_0));
    % Set value of found index to pi
    V(w_0_c) = pi;
    % Find freq index closest to -w_0
    [~, w_0_c] = min(abs(w + w_0));
    % Set value of found index to pi
    V(w_0_c) = pi;
   
    % Calculate approx abs Fourier spectrum based on DFT
    % Scaling is different here, why?
    V_dft = abs(fftshift(fft(v) /  N * 2 * pi));
    
    case 5
    % Combined rect signal -----------------------------
    % Define auxiliary parameters
    a = 1;
    
    % Set signal composed of multiple rect functions
    v = ones(N, 1);                        % Fill signal with ones
    v(t < -T_0) = 0 * v(t < -T_0);     % Set signal to 0 before first rect pulse
    v(t > -T_0) = a * v(t > -T_0);     % Set signal to a from first rect pulse
    v(t > 0) = 2 * v(t > 0);               % Set signal to 2 * a from second rect pulse
    v(t > T_0) = -1/2 * v(t > T_0);    % Set signal to - a from third rect pulse
    v(t > 4*T_0) = 0 * v(t > 4*T_0);  % Set signal to 0 after third rect pulse
    
    
    % Fourier transform according to sample solution
    % Using eps to avoid NaN value in the center
    V = a ./ (1i * w + eps) .* (exp(1i*T_0.*w) + 1 - 3 .* exp(-1i*T_0.*w) + exp(-4i*T_0.*w));

    % Alternative approach to get Fourier transform using rect pulses
    % DSP toolbox only
    % V = 1 * a * T_0 * sinc(w / pi * 0.5 * T_0) .* exp(0.5i * w * T_0);
    % V = V +  2 * a * T_0 * sinc(w / pi * 0.5 * T_0) .* exp(-0.5i * w * T_0);
    % V = V -  3 * a * T_0 * sinc(w / pi * 1.5 * T_0) .* exp(-2.5i * w * T_0);

    % Take absolute value
    V = abs(V);

    V_dft = abs(fftshift(fft(v) /  f_s));
        
end

%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

w_lim = 20;

figure
set(gcf, 'Position', [2000 0 1400 1350]) % Set window size

subplot(3, 1, 1)
plot(t, v, 'LineWidth', 2, 'Color', cmap(1, :))
hold on
grid on
xlabel('$t$ / $\mathrm{s}$', 'FontSize', 20);
ylabel('$v(t)$', 'FontSize', 20);
title('Continuous time signal (discretized)', 'FontSize', 20)

subplot(3, 1, 2)
plot(w, V, 'LineWidth', 2, 'Color', cmap(2, :))
grid on
xlabel('$\omega$ / $\mathrm{s^{-1}}$', 'FontSize', 20);
ylabel('$|V(j\omega)|$', 'FontSize', 20);
title('Fourier spectrum (discretized)', 'FontSize', 20)
xlim([-w_lim, w_lim])

subplot(3, 1, 3)
plot(w, V_dft, 'LineWidth', 2, 'Color', cmap(3, :))
grid on
xlabel('$\omega$ / $\mathrm{s^{-1}}$', 'FontSize', 20);
ylabel('$|V_{\mathrm{dft}}(j\omega)|$', 'FontSize', 20);
title('Fourier spectrum approx. using DFT', 'FontSize', 20)
xlim([-w_lim, w_lim])
