function [DFT_signal, f, magnt, P1, phase] = DFTviafft (Fsignal, Fs)
% DFT computation via FFT algorithm: DFT computation of an the input vector 
% magnitude and phase of the transformed data and plot of the corresponding magnitude and phase spectrum  
%
% @file DFTviafft.m
%
%
% This function applies DFT transformation and performs the following
% operations: 
%  1 - DFT computation
%  2 - Magnitude and phase calculation, of the transformed sequence
%  3 - Plot the magnitude and phase spectrum
%
%
%
%   [DFT_signal, f, P1,P2,phase] = DFTviafft (Fsignal, Fs)
%  
%   INPUTS:
%
%   Fsignal:        Data/Signal to be processed. Signals should be filtered
%                   previously if that step is required.  
%                   * A Butterworth filter is provided on @file Butt_highpass.m
%
%   Fs:             Sampling Frequency    
%
%  OUTPUTS:
%
%  DFT_signal:      DFT of the input vector 
%
%  f:               Frequency vector  
%
%  magnt:           Magnitude of the transformed sequence
%
%  P1:              Single-Sided amplitude spectrum of the transformed sequence
%
%  phase:           Phase of the transformed sequence
%
%
% Beatriz Jacinto Barros (beatriz.j.barros@inesctec.pt)
% BRAIN, INESC-TEC Porto, October 2022
% 

%% Calculate the parameters of the signal
T = 1/Fs;                                           % Sampling period
L = length(Fsignal)                                 % Length of signal  %signal: time-domain signal
t = (0:L-1)*T;                                      % Time vector    

%% Compute Discrete Fourier transform (DFT) of X using a fast Fourier transform (FFT) algorithm
DFT_signal= fft(Fsignal);                           % Fourier transform of the signal
df=Fs/L;                                            % Frequency resolution
sampleIndex = -L/2:L/2-1;                           % Ordered index for FFT plot
f=sampleIndex*df;                                   % X-axis index converted to ordered frequencies

%% Compute magnitude and phase of the transformed sequence

magnt=abs(DFT_signal);                                      % Magnitude

% Decrease round-off error when computing the phase by setting small-magnitude transform values to zero
DFT_signal(magnt<1e-6) = 0;

% Phase computation with phase unwrapping algorithm 
% compensate  extrinsic’ discontinuities on the phase spectrum
phase = unwrap(angle(DFT_signal));                               % Phase

%% Plot phase (Degrees)
figure(1)
plot(f,phase);                                      % Plot Phase vs Frequencies
ylabel('Phase [Degrees]'), grid on
xlabel('Frequency [Hertz]') 

%% Two and Single-Side Spectrum
P2 = abs(DFT_signal/L);                             % Two-sided spectrum
P1 = P2(1:L/2+1);                                   % Single-sided spectrum
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;                                 % Frequency Domain 

%% Plot amplitude spectrum
figure(2)
plot(f, P1)                                         % Plot Magnitudes vs Frequencies
title("Single-Sided Amplitude Spectrum of X(t)")
xlabel("f (Hz)")
ylabel("|P1(f)|")





